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Preface 


Mathematical modelling is the process of replacing problems from outside mathemat- 
ics by mathematical problems. The subsequent mathematical treatment of this model 
by theoretical and/or numerical procedures works as follows: 


1. Transition from the outer-mathematics phenomenon to a mathematical descrip- 
tion, what at the same time leads to a translation of problems formulated in terms 
of the original problem into mathematical problems. 


This task forces the scientist or engineer who intends to use mathematical tools 
to 


— cooperate with experts working in the field the original problem belongs to. 
Thus, he has to learn the language of these experts and to understand their 
way of thinking (teamwork). 


— create or to accept an idealized description of the original phenomena, i.e. an 
a-priori-neglection of properties of the original problems which are expected 
to be of no great relevance with respect to the questions under consideration. 
These simplifications are useful in order to reduce the degree of complexity 
of the model as well as of its mathematical treatment. 


— identify structures within the idealized problem and to replace these struc- 
tures by suitable mathematical structures. 


2. Treatment of the mathematical substitute. 
This task normally requires 
— independent activity of the theoretically working person. 
— treatment of the problem by tools of mathematical theory. 


— the solution of the particular mathematical problems occuring by the use of 
these theoretical tools, i.e differential equations or integral equations, op- 
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8 Preface 


timal control problems or systems of algebraic equations etc. have to be 
solved. Often numerical procedures are the only way to do this and to answer 
the particular questions under consideration at least approximately. The er- 
ror of the approximate solution compared with the unknown so-called exact 
solution does normally not really affect the answer to the original problem, 
provided that the numerical method as well as the numerical tools are of suf- 
ficiently high accuracy. In this context, it should be realized that the question 
for an exact solution does not make sense because of the idelizations men- 
tioned above and since the initial data presented with the original problem 
normally originate from statistics or from experimental measurements. 


3. Retranslation of the results. 


The qualitative and quantitative statements received from the mathematical mod- 
el need now to be retranslated into the language by which the original problem 
was formulated. This means that the results have now to be interpreted with 
respect to their real-world-meaning. This process again requires teamwork with 
the experts from the beginning. 


4. Model check up. 


After retranslation, the results have to be checked with respect to their relevance 
and accuracy, e.g. by experimental measurements. This work has to be done 
by the experts of the original problem. If the mathematical results coincide 
sufficiently well with the results of experiments stimulated by the theoretical 
forecasts, the mathematical part is over and a new tool to help the physicists, 
engineers etc. in similar situations is born. 


Otherwise — if no trivial logical or computational errors can be found — the model 
has to be revised. In this situation, the gap between mathematical and real results 
can only originate from too extensive idealizations in the modelling process. 


The development of mathematical models does not only stimulate new experiments 
and does not only lead to constructive-prognostic — hence technical — tools for physi- 
cists or engineers but is also important from the point of view of the theory of cog- 
nition: It allows to understand connections between different elements out of the un- 
structured set of observations or — in other words — to create theories. 


Fields of applications of mathematical descriptions are well known since centuries in 
physics, engineering, also in music etc. In modern biology, medicine, philology and 
economy mathematical models are used, too, and this even holds for certain fields of 
arts like oriental ornaments. 
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This book presents an introduction into models of fluid mechanics, leads to important 
properties of fluid flows which can theoretically be derived from the models and shows 
some basic ideas for the construction of effective numerical procedures. Hence, all 
aspects of theoretical fluid dynamics are addressed, namely modelling, mathematical 
theory and numerical methods. 


We do not expect the reader to be familiar with a lot of experimental experiences. 
The knowledge of some fundamental principles of physics like conservation of mass, 
conservation of energy etc. is sufficient. The most important idealization is — contrary 
to the molecular structure of materials — the assumption of fluid-continua. 


Concerning mathematics, it can help to understand the text more easily if the reader 
is acquainted with some basic elements of 


Linear Algebra 


Calculus 


Partial Differential Equations 


Numerical Analysis 


Theory of Complex Functions 


Functional Analysis. 


Functional Analysis does only play a role in the somewhat general theory of dis- 
cretization algorithms in chapter 6. In this chapter, also the question of the existence 
of weak entropy solutions of the problems under consideration is discussed. Physicsts 
and engineers are normally not so much interested in the treatment of this problem. 
Nevertheless, it had to be included into this presentation in order not to leave this 
question open. The non-existence of a solution shows immediately that a model does 
not fit the reality if there is a measurable course of physical events. Existence theo- 
rems are therefore important not only from the point of view of mathematicians. But, 
of course, readers who are not acquainted with some functional analytic terminology 
may skip this chapter. 


With respect to models and their theoretical treatment as well as with respect to nu- 
merical procedures occuring in sections 4.1, 5.3, 6.4 and chapter 7, a brief introduction 
into mathematical fluid mechanics like this book can only present basic facts. But the 
author hopes that this overview will make young scientists interested in this field and 
that it can help people working in institutes and industries to become more familiar 
with some fundamental mathematical aspects. 
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Finally, I wish to thank several colleagues for suggestions, particularly Thomas Sonar, 
who contributed to chapter 7 when we organized a joint course for graduate students!, 
and Dr. Michael Breuss, who read the manuscript carefully. Last but not least I thank 
the publishers, especially Dr. Alexander Grossmann, for their encouragement. 


Rainer Ansorge 


Hamburg, September 2002 


‘Parts of chapters 1 and 5 are translations from parts of sections 25.1, 25.2, 29.9 of: Ansorge and 
Oberle: Mathematik fiir Ingenieure, vol. 2, 2nd ed., Berlin, Wiley-VCH 2000. 
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1 Ideal Fluids 


1.1 Modelling by Euler’s Equations 


Physical laws are mainly derived from conservation principles, such as conservation 
of mass, conservation of momentum, conservation of energy. 


Let us consider a fluid (gas or liquid) in motion, i.e. the flow of a fluid’, let 


ui (2, Y, Z, t) 
u(x, y, z, t) = u(x, Y, Z, t) 
u3(z, Y, Z, t) 


be the velocity and denote by p = p(x,y,z,t) the density of this fluid at the point 
x = (x,y,z) and at the time instant t?. 


Let us take out of the fluid at a particular instant t an arbitrary portion of volume W (t) 
with surface OW (t). The particles of the fluid now move, and assume W(t + h) to 
be the volume at the instant t + h formed by the same particles which had formed 
W (t) at the time t. 


Moreover, let y = y(x, y, z,t) be one of the functions describing a particular state of 

the fluid at the time ¢ at the point x such as mass per unit volume (= density), interior 

energy per volume, momentum per volume etc. Hence, f yd(zx,y,z) gives the 
W(t) 

full amount of mass or interior energy, momentum etc. of the volume W(t) under 

consideration. 


‘Blows of other materials are eventually included, too, e.g. flow of cars on highways, provided that 
the density of cars or particles is sufficiently high. 
Bold letters normally mean vectors or matrices. 
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14 1 Ideal Fluids 


We ask for the change of f yd(«,y,z) with respect to time, i.e. for 


W(t) 
d 
dt f p (x,y, z,t) d(£,y, z). (1.1) 
W(t) 
We have 
z | ovatae) 
dt p T, Y, Z, T, Y, z 
W(t) 
. 1 E 
E J p(Ņ, t+ h)d(y1, y2, Y3) — f p(x, y,z,t)d(x, y, z) , 


W(t+h) W(t) 
where the change from W(t) to W(t + h) is obviously given by the mapping 
y= x2+h-u(a,t) + o(h) 


(y = (yi, y2,93)-). 


The error term o(h) also depends on æ but keeps the property lim 70(h) = 0 if 


differentiated with respect to space provided that these spatial derivatives are bounded. 


The transformation of the integral taken over the volume W (t + h) to an integral over 
W (t) by substitution requires the integrand to be multiplied by the determinant of this 
mapping, i.e. by 


(1 + h Oru) hOyur hozuy 
h zuz (1+ hdyu2) h ðzu2 + o(h) 
h ðru3 h Oyug (1 + h ðzu3) 


= 1 + h (Ozu, + Oyu2 + O,u3) + o(h) 


= 1 + h- divu(z,t) + olh). 
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Taylor expansion of y(y,t + h) around (x,t) therefore leads to 


d 
dt I p(x, y, Z,t) d(x, y, z) 


W(t) 
= | {ap+edvas(uve)} dzy). 02) 
W (t) 


Herewith, Vv denotes the gradient of a scalar function v, and (- , -) means the stan- 
dard scalar product of two vectors out of R3. 


The product rule of differentiation gives 
pdivu+(u,Vy) = div(y-u), 
such that (1.2) leads to the so-called Reynolds’ transport theorem? 


d 
i i g(x,y, 2,t)d(x,y,z) = f {hp + div(pu)} d(z,y,z). (1.3) 


Wit) W(t) 


As already mentioned, the dynamics of fluids can be described directly by conser- 
vation principles and — as far as gases are concerned — by an additional equation of 
state. 


1. Conservation of mass: If there are no sources or losses of fluid within the sub- 
domain of the flow under consideration, mass remains constant. 


Because W(t) and W(t + h) consist of the same particles, they are of the same mass. 
The mass of W(t) is givenby f p(x,y,z,t) d(x, y, z), and therefore 
W(t) 


£ f p(x, y, z,t) d(x,y, 2) = 0 
W(t) 


must hold. Taking (1.3) into account (particularly for y = p), this leads to the 
requirement 


{2p + div(pu)} d(z,y,z) = 0. 
W(t) 





3Osborne Reynolds (1842 - 1912); Manchester 
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Since this has to hold for arbitrary W(t), the integrand has to vanish: 


Op + div (pu) = Ol (1.4) 


This equation is called the continuity equation. 


2. Conservation of momentum: Another conservation principle concerns the mo- 
mentum of a mass which is defined as 


mass x velocity. 


Thus, 
f pud(z,y,z) 
W(t) 


gives the momentum of the mass at time t of the volume W(t) and 


describes the density of momentum. 
The principle of conservation of momentum, i.e. Newton’s second law 


force = mass x acceleration, 


now states that the change of momentum with respect to time equals the sum of all 
exterior forces acting on the mass of W (t). 


In order to describe these exterior forces, we take into account that there exists a 
certain pressure p(x,t) at every point a of the fluid and at every instant t. If n 
is considered to be the unit vector normal on the surface OW(t) of W(t) and of 
outward direction, the fluid outside of W(t) leads to a force acting on W(t) which 
is given by 


— / pndo (do = area element of OW(t)) . 
awit) 


Besides the normal forces per unit of the surface area generated by the pressure, there 
exist also tangential forces which act on the surface by friction of the exterior particles 
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along the surface. Though this so-called viscosity of a fluid leads to a lot of remark- 
able phenomena, we are going to neglect this property at the first step: Instead of real 
fluids or viscous fluids we restrict ourselves in this chapter to so-called ideal fluids 
or inviscid fluids. This restriction to ideal fluids, particularly to ideal gases, is one of 
the idealizations mentioned in the preface. 


But there are not only exterior forces per unit of the surface area but also exterior 
forces per unit volume, e.g. the weight. 


Let us denote these forces per unit volume by k such that Newton’s second law leads 
to 

d 

dt qd(x, y, z) = k (x,y, 2, t) d(x,y, z) g p:-ndo. 


W(t) W(t) OW (t) 


Thus, by Gauss’ divergence theorem, we find 


pn, do J O,p d(x, y, z) 
aw (t) W(t) 
/ pies f pnzdo | _ f Oyp a(x, y, z) 
IWA) aw (t) (t) 
png do / dzp d(x,y, 2) 
aw (t) W(t) 


= f Vpd(z,y,z). 


W(t) 
Together with (1.3), 
div (qı u) 
i ðq+ | div(qau) | -—k+Vp> da,y,z) = 0 
wit) div (q3 u) 
follows. 


Again, this has to be valid for every arbitrarily chosen volume W(t). If, moreover, 
div (qu) = (u,Vq) + divu-G 
is taken into account, 


dq + (u,V)qt+divu-q+ Vp = k, 
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1 1 
Cg dy Nal ogy (ža) a+ vp = k (1.5) 





has to be fulfilled. 


The number of equations represented by (1.5) equals the spatial dimension of the flow, 
i.e. the number of components of q or u. 


By means of the continuity equation, (1.5) can be reformulated as 





jut (u,V)u+iVp = k, (1.6) 
where the force k per unit volume has been replaced by the force k = 7 k per unit 
mass. 


Equation (1.4) together with (1.5) or (1.6) are called Euler’s equations’. 


P, E, qi, 92, 43 


are sometimes called conservative variables whereas p, €, u1, u2, ug are the prim- 
itive variables. Herewith, 


P42 llall? 
E := pra = 
pet 5 Ilall pE + Zp 


gives the total energy per unit volume where £ stands for the interior energy per 
unit mass, e.g. heat per unit mass. 


P 2 
lul 


obviously introduces the kinetic energy per unit volume. 


(u, V ) u is called the convection term. 


Remark: D 
Terms of the form ðw + ( w, V ) w are in the literature often abbreviated by aoe 


and are called material time derivatives of the vector valued function w. 


4Leonhard Euler (1707 - 1783); Basel, Berlin, St. Peterburg 
5]| - || : 2-norm. 
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Next we consider the First Law of Thermodynamics, namely the 


3. Conservation of energy. 


The change per time unit of the total energy E of the mass of a moving fluid volume 
equals the work which is done per time unit against the exterior forces. 


This means in the case of ideal fluids 


d 
dt f E(x, y, z,t) d(x, y, z) 


W(t) 


The relation 


(pnu)do = | (pun) do = f div(puyd(e.y.2) 


ƏW (t) aw (t) W(t) 
follows from the Gauss’ divergence theorem, so that (1.3) leads to 


J {aE + div (E - u) + div (p- u) — (pk, u) }d(z,y,2) =0 VWi(t). 
W(t) 


If k can be neglected because of small gas weight or if k is the weight of the fluid 
per unit mass’ and u the velocity of a flow parallel to the earth surface, we get 


AE + div (= a) = o0]. AT) 





work 
time 





ĉwork = force x length = pressure x area x length —> = force x velocity = pressure x 
area x velocity. 


Tie. ||k|| = g with the earth acceleration g. 
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Explicitly written and neglecting k , the equations (1.4), (1,5) and (1.7) read as 




















dtp + Ord + 0yq2 + 0.93 =0 
1 1 1 

Ong + Ox G ag +p) + Oy G ng) + Oz G ng) =0 
p p p 
1 1 1 

Ong2 + Ox (- on) + Oy (- q2q2 +p) +0, (- 2a) =0 
p p p 
1 1 1 

193 + Oz C on) + Oy C wats) + 0z G 9393 +p) =0 
p p p 
E+ E+ E+ 

B+ De ( 2a) +0, ( =? an) +a ( =P as) z3 


Hence, we are concerned with a system of equations to be used in order to determine 
the functions p, q1, q2, q3, E. But we have to realize that there is an additional 
seeked function, namely the pressure p. In case of constant density pÈ, i.e. Op = 
0 Y(z,y,z), only four conservation variables have to be determined such that the 
five equations (1.8) are sufficient. Otherwise, particularly on event of gas flow, a sixth 
equation is needed, namely an equation of state. State variables of a gas are 


T = temperature, p = pressure, p = density, V = volume, 


€ = energy/mass, S = entropy/mass , 


and a theorem of thermodynamics says that each of these state variables can uniquely 
be expressed in terms of two of the other state variables. Such relations between three 
state variables are called equations of state. Thus, p can be expressed by p and e€ 
(hence by p, E, q) ; for inviscid so-called y-gases the relation is given by 


2 
q 
p= 0- Dee = (v-1) (B- L) 19) 
p 
with y = const > 1, such that only the functions p, q, Æ have to be determined. 
Herewith, ~y is the ratio = of the specific heats. If air is concerned, we have 
y = 1.4. 


By means of vector valued functions, also systems of differential equations can be 
described by a single differential equation. In this way and taking (1.9) into account 


The case of constant density is not necessarily identical with the case of an incompressible flow 
defined by divu = 0 because the continuity equation (1.4) is in this situation already fulfilled if the 
pair (p , u) shows the property Ojo + lu , Vp) = 0, which does not necessarily imply p = const. 
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as an additional equation, (1.8) can be written as 





OV + Or fi(V) + 3y foV) + Oz fV) = 0 (1.10) 
with 
T 
V = (p, qi, q2, q3, E) 
and with 
qı q2 q3 
1 1 1 
— qıqı +P — Qa — Bq 
pP P P 
1 1 1 
KSI Gen- prys poe Vi |) ee 
1 1 1 
— 3% — 93q2 — 9393 + p 
P P p 
E+p E+p E+p 
— q —_ 42 q3 
P P P 


The functions f;(V) are called fluxes. 


If J fi, J fo, J fz} are the Jacobians, respectively, (1.10) reads as 


Because of their particular meaning in physics, systems of differential equations of 
type (1.10) are called systems of conservation laws, even if they do not originate 
from physical aspects. Obviously, (1.11) is a quasilinear systems of partial differential 
equations of first order. 


Such systems have to be completed by initial conditions 
V(x,0) = Vo(z) , (1.12) 
where Vo(a) is a prescribed initial state, and by boundary conditions. 


If the flow does not depend on time, it is called a stationary or steady-state flow and 
initial conditions do not occur. Otherwise, the flow is called an instationary one. 


22 1 Ideal Fluids 


As far as ideal fluids are concerned, it can be assumed that the fluid flow is tangential 
along the surface of a solid body? fixed in space or along the bank of a river etc. In 
this case, 


(u,n) =0 (1.13) 


is one of the boundary conditions where n are the outward normal unit vectors along 
the surface of the body. 


Often there occur symmetries with respect to space such that the number of unknowns 
can be reduced, e.g. in the case of rotational symmetry combined with polar coordi- 
nates. If there remains only one spatial coordinate, the problem is called one dimen- 
sional!°. If rectangular space variables are used and if x is the only remaining one, 
we end up with the system 


aV + Jf(V)-.V = 0 (1.14) 
with 
q 
P 12 
V = = pd 
l ; f(V) ae 
—— 4 
p 


and with the equation of state 


p = (y-!) (e-£). 


2p 
Hence, 
0 1 0 
sag q 
Jf(V) = oe ge PEN y-1 
3 2 
q Eq 3(y-1) q q 
GaP e T a 
p p 2 P p 


°e.g. of a wing 
‘Two or three dimensional problems are defined analogously. 
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As can easily be verified, A; = 1 is a solution of the characteristic equation of 
p 
J f(V) namely of 


2 E 
-43i S-a nE) 


2p? 
Y y q q 
+(S-341) pres eee 0 


The other two eigenvalues of J f(V) are then the roots of 


2 8 E 
i MS Ia) Se 




















2p? p 
2 2 
q q E q 
à = 2a tya- =- aea) a 
p pP p 2p? 
= t4 ho- = - vy D Gaa 
p p 2p 
—1 2 
a ie O |e- E] = £4 E 
p p 2p p p 


Thus, these eigenvalues are real and different from eachother (p > 0). 


Definition: 
If all the eigenvalues of a matrix A(x,t, V) are real and if the matix can be dia- 
gonalized, a system of equations 


AV + A(z,t,V)0,V = 0 


is called of hyperbolic type at (x,t,V). If the eigenvalues are real and different 
from eachother so that A can certainly be diagonalized, the system is called strictly 
hyperbolic. 


So, obviously, (1.14) is strictly hyperbolic for all (x,t, V) under consideration. 


By the way, because of r= u, the velocity of the flow, and because , / IP describes 
P P 


the local sound velocity ¢ (cf. (1.67)), the eigenvalues are 


à = u, Ag = ute, 43 = u-é, 
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and have equal signs in case of supersonic flow whereas the subsonic flow is charac- 
terized by the fact that one eigenvalue shows a different sign than the others. 


Definition: 

Let u(x, to) be the velocity field of a flow at the instant to and let xo be an arbitrary 
point of the particular subset of R? which is occupied by the fluid at this particular 
instant. If at this moment the system 


[æ"(s) , u(w(s),to)] = 07 
of ordinary differential equations with the initial condition 
(0). = £o 
has a unique solution æ = a(s) for each of these points æo, each of the curves 


x = g(s) is called a streamline at the instant to. 


Here, the streamlines may be parametrised by the arc length s with s = 0 at the 
point £o, such that 
(a a’) =z 


Thus, the set of streamlines at an instant tọ shows a snapshot of the flow at this 
particular instant. It does not necessarily describe the trajectories along which the 
fluid particles move when time goes on. Only if the flow is a stationary one, i.e. if 
u, p, P, k are independent of time, streamlines and trajectories coincide: A 
particle moves along a fixed streamline in the course of time. 


1.2 Characteristics and Singularities 


As an introductory example for a more general investigation of conservation laws, let 
us consider the scalar Burgers’ equation! without exterior forces 


1 
Ow + Op (507) =0. (x,t) EQ, iexeR,t>0, (1.16) 


which is often studied as a model problem from a theoretical point of view and also 
as a test problem for numerical procedures. 


"where [- , -] is the vector product in R? 
12], Burgers: Nederl. Akad. van Wetenschappen 43 (1940), 2-12. 
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Here, the flux is given by f(v) = 507. If 


x 

volz) =1— 5 (1.16a) 
is chosen as a particular example of an initial condition, the unique and smooth so- 
lution!’ turns out to be 


2-2 


v(x, t) = a 





but this solution does only exist locally, namely for 0 < t < 2; for increasing time it 
runs into a singularity at t = 2. 


As a matter of fact, classical existence and uniqueness theorems for quasilinear par- 
tial differential equations of first order with smooth coefficients ensure the unique 
existence of a classical smooth solution only in a certain neighbourhood of the initial 
manifold, provided that also the initial condition is sufficiently smooth. 


The occurence of discontinuities or singularities does not depend on the smoothness 
of the fluxes: Assume v(x,t) to be a smooth solution of the problem 


Ou + ðf) =0 for rER,t>0 


v(x, 0) = volz) 


in a certain neighbourhood immediately above the x-axis. Obviously, this solution is 
constant along the straight line 


a(t) = zo +t f'(vo(xo)) (1.17) 


that crosses the x-axis at xo where xo € R is chosen arbitrarily. Its value along this 
line is therefore u(x(t),t) = vo(xo). This can easily be verified by 





These straight lines (1.17), each of which belonging to a particular xo, are called the 
characteristics of the given conservation law. 





'3; e. v is continuously differentiable 
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th 


u = u (x1) 











= 
X0 X1 X 


Figure 1: Formation of discontinuities 


In case of xo < 21, but e.g. 0 < f'(vo(£1)) < f'(volzo))!t, the characteristic 
through (xo,0) intersects the characteristic through (x1,0) at an instant tı > 0, so 
that at the point of intersection the solution v has to have the value vo(zo) as well as 
the value vo(zı) Æ vo(xo). Therefore, a discontinuity will occur at the instant tı or 
even earlier. With respect to fluid dynamics, particularly shocks — discussed later on 
— belong to the various types of discontinuities. 


Applied to Burgers’ equation (1.16) with the initial function (1.16a), the characteristic 
through a point xo of the x-axis is given by 


T — LH 





t=2 : 
2 — To 


such that the discontinuity we had found at t = 2 can also immediately be understood 
by means of Figure 2. 


If systems of conservation laws are considered instead of the scalar case, i.e. 
V(x, t) ER”, mEN,Y(z, t) EQ, 


and if only one spatial variable x occurs, a system of characteristics is defined as 
follows: 


Definition: 
If V(a,t) is a solution of (1.10) in the case of only one space variable x , the 
one-parameter set 


Ta = Talt, x) (x € R : set parameter) 





l4hence also vo(zo) Æ vo(x1) 
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x 


Figure 2: Discontinuity of the solution of the Burgers equation 


of real curves defined by the ordinary differential equation 


t = A(V(x;2)) (1.18) 
for every fixed i (i =1,...,m) is called the set of i-characteristics of the partic- 
ular system that belongs to V. Here, \;(V) (i= 1,---,m) are the eigenvalues 


of the Jacobian J f(V). 


Obviously, this definition coincides in case of m = 1 with the previously presented 
definition of characteristics. 


Let us finally and exemplarily study the situation of a system of conservation laws if 
the system is linear with constant coefficients: 


OAV + AV =0 (1.19) 


with a constant Jacobian (m, m)-matrix A. Moreover, let us assume the system to be 
strictly hyperbolic. 
From (1.18), the characteristics turn out to be the set of straight lines given by 


zalt) = At+2@(0), (=1,---m), 


and independent of V. Hence, for every fixed 2, the characteristics belonging to this 
set are parallel. 


If S = (81, 82,--+ , 8m) is the matrix whose columns consist of the eigenvectors of 
the Jacobian and if A denotes the diagonal matrix consisting of the eigenvalues of A, 


A = SAS"! 


follows. 
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If new variables V are introduced by 
v=s'v , 
the sytem gets the form 


AV +A0, =0, Vo =S Vo. (1.21) 


This is a decoupled system: 
O10; + ri Ox0; = 0 5 ;(x, 0) = [S~'Vo(2)], = Vig (2) 5 (i = 1, ps ,m) fs 


Each of the equations of this system is an independent scalar equation called an ad- 
vection equation. The solution is 


ilz, t) = on Oa a Ài t) (i _ 1, see ,m) . (1.22) 


Hence, the state at the instant t moves with velocity A; into positive or negative x— 
direction according to the sign of A;, respectively. This is called a wave propagation 
where the propagation velocity is described by A;. 


Obviously, 
V = (sı sie Sp) : = 018, + to82 +... + Um Sm. 


Each of the vector valued functions 
Vw (a, t) := CEL) -> si 
solves the system of differential equations because of 


Va + SAST IVa = Odi si + Ordi Ai si 


(dôi + Ai Oz 0;) g = O 


Va (4 =1,---,m) is often called the solution belonging to the i-th set of char- 
acteristics and to the given initial value ĉ;(x). Obviously, the vector functions 
Va (i=1,:-:,m) are linearly independent. 

Remark: 


The fact that a sufficiently smooth solution of a nonlinear initial value problem does 
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often only exist in the neighbourhood of the initial manifold —though the correspon- 
dent real world process shows a global existence— can’t be accepted by a physicist 
or engineer. He insists the mathematician to present a global solution. This forces 
the mathematicians to create a more general definition of a solution such that the real 
world situation can be described in a satisfactory way. Suitable definitions of weak 
solutions will be presented in chapter 2. 


Remark: 
As far as the scalar linear problem 


Ow + 0,(av) =0 , a=const , (1.23) 
is concerned, v often describes a concentration, and the flux then is simply given by 
flv) = an. 


Many physical processes include a further flux of the particular form 
—cdO,v (e>0) , 


proportional to the concentration descent. The transport phenomenon is then mathe- 
matically modelled by 


Ow + Oy (av—E0,v) = 0, ie. Ov + ad,v = EOzzv. (1.24) 


Because of the diffusion term on the right hand side, this equation is of parabolic 
type. Parabolic equations make their solutions smoother and smoother when time t 
increases. Therefore, shocks will be smeared as soon as diffusion occurs. 


This effect of parabolic equations can easily be demonstrated by examples like the 
following one: The sum 


OO . 
v(t) = ef) epi ioa 
(x,t) Lear 
converges uniformly for t > 0 because of the factors e~” *t_ This also holds after 
multiple termwise differentiations of this sum with respect to t as well as with respect 
to x so that v is a sufficiently smooth function for all x € R , t > 0. It solves 
(1.24)! in the classical understanding, but leads for t = 0 to the function 





Sin case of £ = 1 
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where w(x) represents the Fourier expansion of the 27-periodic discontinuous func- 
tion 





(j= 5 for -7<au<7 
Cet 0 for v= ET 


The curve described by w(x) is sometimes called a saw blade curve. 
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Let us try to investigate the forces acting on solid bodies!® dipped into a fluid flow at 
a fixed position. For convenience, we restrict ourselves to stationary flows of inviscid 
fluids. Moreover, it will be assumed that the magnitudes of the velocities are such 
that the density can be regarded as a costant. Thus, the flow is incompressible and the 
partial derivatives with respect to t occuring in (1.8) vanish. The first equation in (1.8) 
if written in primitive variables, i.e. the continuity equation (1.4), then reduces to 


divu = 0, (1.25) 


T17 


where u = (u(x, y, z), ua(x,y, z), uz3(z,y,z)) again stands for the velocity 


vector of the flow at the space position x = (x,y,z) . 


The second, third, and fourth equation of (1.8) formulated by means of primitive 
variables could be written as the vector valued Euler equation (1.6) and lead in case 
of a stationary flow to 


ues a = k. (1.26) 
p 


Let us —moreover— assume that the flow is irrotational; this means that the circula- 
tion 


LZ. = f ula) de (1.27) 


vanishes for every closed contour C within every simply connected subdomain of the 
flow area. Then w can be derived from a potential @ i.e. there is a scalar function 
@ so that 


u=Vo, (1.28) 


‘e.g. on a wing of an aircraft 
For the time being, z means the third space variable; it will later denote the complex variable 
x + iy, but confusion will certainly not occur. 
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and 
curl u = 0 (1.29) 


holds within this area. 


By the way, the vector curl u is often called the vorticity vector or angular velocity 
vector, and a trajectory of a field of vorticity vectors is called a vortex line. 


Because of (1.28), a flow of this type is called a potential flow, and ¢ is the so-called 
velocity potential. 


The fifth equation in (1.8) can be omitted as far as the knowledge of the energy density 
E is of no interest. 


From (1.29), the relation (u,V)w=4V/((lull?)—[u, curl u] leads together 
with (1.26) to 


1 1 A 
v ($lel? + +p) Lw (1.30) 
p 


Formula (1.30) shows that a fluid flow of the particular type under consideration, 
namely an approximately stationary, inviscid, incompressible, and irrotational flow 
can only exist if the exterior forces are the gradient of a scalar function. In other 
words: these forces have to be conservative, i.e. there has to exist a potential Q with 
k= -V Q, such that (1.30) leads to 


1 1 
v (5l? +2p+Q) = 0. 
p 
i.e. to 
Neca | 
z lll + a +Q = const. (1.31) 


(1.31) is called the Bernoulli equation! and is nothing else than the energy conser- 
vation law for this particular type of flow. p is called the static pressure, whereas 
the term £ \|20||? i.e. the kinetic energy per volume, is often called the dynamic 
pressure. 


As far as an incrompressible flow in a circular pipe is concerned, the velocity will 
necessarily increase as soon as the diameter of the pipe decreases, and — because of 





'8Daniel Bernoulli (1700 - 1782); St. Peterburg, Basle 
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(1.31) — this will lead to a decreasing pressure within the narrow part of the pipe. This 
phenomenon is called the hydrodynamic paradoxon. Applications are carburettors, 
jet streams etc. 


Remark: 

A necessary condition for the existence of the irrotational flow was the conservative 
character of the exterior forces. Let us now assume that —vice versa— these forces 
are conservative, i.e. k= -V Q. Morover, let us admit that the flow is even 
compressible in the sense of a particular dependence p = p(p), called a barotropic 
flow. Integration of (1.26) along a streamline from a constant point Po to a variable 
point P then leads to 


P. 
READE E E et vo} din 


2 
Po 
with 
d 
O -= | £ ; 
p(p) 
hence 
VƏ = -Vp , 
p 
and with 
d d ( length) 
s = — uds s = arc length) . 
lull 


Because P was arbitrary and because of (|curlu, u], u) = 0, this result leads to 
the generalized Bernoulli equation 


1 
3 lul? + © +Q = const, 


which in the case of constant density, i.e. © = 2 seems to coincide completely with 


(1.31). But one has to take into account that the constant at the righthand side may 
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now change from streamline to streamline. Additionally, we find 


d d d da 
a = gf ouf Glue pe 
C 


C 
du 
= f Zag udu 
C C 
du 
Z f Zag (lap Wal de 
dt 
C C 
1 


= -$ V(0+Q) de =o 
C 


This is Kelvin’s Theorem!”, which says that the circulation along a closed curve in 
an inviscid barotropic flow does not change in time. 


If irrotational flows of arbitrary type of incompressibility are under consideration, 
(1.28) leads together with (1.25) to 


Ag=0 (1.32) 


whereas the continuity equation (1.4) reads for compressible fluids in the case of sta- 
tionary flows as 


div(pu) =0. 
Let us additionally assume the flow to be barotropic. Then, because of 
; f dp 
div(p u) = p div u + — (Vp , u) 
dp 
and because of 


ga (cf. (1.68)) , 
dp 


(1.26) leads to 














Lord Kelvin of Largs (1824 - 1907); Glasgow 
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as far as the exterior forces vanish. In this situation, (1.33) generalizes (1.32). 
Obviously, (1.33) is a quasilinear partial differential equation of second order for œ, 
certainly elliptic if 


M := leli <1, (1.34) 
C 


i.e. in areas of subsonic flow. 


Definition: M is called the Mach number”. 


If particularly a constant flow into x-direction is under consideration which is only 
disturbed in the neighbourhood of a slim airfoil?! with a small angle of attack, prod- 
ucts of the values u; (i = 2, 3) with each other and with uw; can be neglected 
compared with 1. This leads to ||u||? = u1? and to a shortened version of (1.33), 
namely to 


U1 


A 
(1 7 (=) ) na + Oy + O26 =0 . (1.35) 


This equation is hyperbolic in areas of supersonic flow (W > 1), and also the full 
equation (1.33) is of hyperbolic type in this case. 





a 
eee 
Figure 3: Flow around a wing if the angle of attack is small 


Let us extend our idealizations by the assumption that the flow under consideration is 
a two-dimensional plane flow. This means that one of the components of the veloc- 
ity vector u in a rectangular coordinate system, e.g. the component in z-direction, 


Ernst Mach (1838 - 1916); Graz, Praha, Vienna 
>! cross-section of a wing or another rigid body in a plane parallel to the direction of the flow 
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vanishes for all (x, y) € R? and for all t > 0: 
u3 = 0. (1.36) 


In case of a two-dimensional supersonic flow along a slim airfoile for which (1.35) 
holds, Mach’s angle 8, determined from 


1 
; E 
| sin 6| M S 
describes the angle between the characteristics of the wave equation 
(1 = M?) dre + yy =0 


and the flow direction given by the direction of the x-axis. The set of all these char- 
acteristics is called Mach’s net. This net plays an important role as far as so-called 
methods of characteristics are used in order to establish efficient numerical proce- 
dures. 








Figure 4: Mach’s net in the case of a linearized supersonic flow (M = const > 1) 


The two-dimensional case of rotational symmetry, e.g. flow along a projectile, leads 
correspondingly to Mach’s cone. 


In the two-dimensional plane situation, (1.29) reads as 
—0,uU2 
curl u = zu = 0: 


O,U2 = Oyu 
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xT 











Figure 5: Linearized supersonic flow around a slim rotatory cone 


Thus, u and ug are independent of the third spatial variable z : 


ulz, y) 
u = u2(2, y) 
0 
with 
Or U2 = Oyu =" 0 (1.37) 


After introduction of the vector 
—u2(z, y) 
v := u(x, y) (1.38) 
0 


and by means of the continuity equation (1.25), we receive 


0 0 
curl v = 0 = 0 = O, 
Ozu + yuz div u 


such that also v can be derived from a potential in simply connected parts of the fluid 
area. In other words, there is a scalar function Y = w(x, y) with 


—uz = Oz, UL yy) (1.39) 


yw is called the stream function. 
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(1.28) leads to 
Ore = Oy, Oyo = —Oz~p. (1.40) 


Obviously, (1.40) can be interpreted as the system of Cauchy—Riemann equations 
of the complex function 


Q(z) := P(az,y) + iyp(z, y) (1.41) 
depending on the complex variable z = x + iy. 


Q is called the complex velocity potential of the plane potential flow under consid- 
eration. 


We are going to assume that the first partial derivatives of the functions @ and w 
are continuous such that Q turns out to be a holomorphic function, and we intent to 
study the forces acting on rigid bodies when dipped into such a fluid flow. 


Because we reduced the reality to a plane flow, the rigid body is assumed to be very 
long with respect to the direction of the third spatial variable; more precisely, it has to 
be of infinite length from the point of view of mathematics. Therefore, the flow around 
the contour of the cross section of the body within the (x, y)-plane, e.g. around the 
contour of an airfoile of a long wing, is of interest. 


We stated that the circulation Z = $ udæ vanishes as long as the closed contour 


C 
C is the boundary of a simply connected domain within the fluid area in the case of a 
potential flow. Because of our plane model, C represents a simply closed contour in 
the (x, y)-plane, i.e. in R°. 


But now, a rigid body with its boundary IĮ is dipped into the fluid, e.g. a wing. If its 
airfoil is part of the area surrounded by C', this interior domain of C is no longer a 
simply connected domain of the fluid area. Hence, the circulation Z around the airfoil 
does not necessarily vanish but turns out to fulfill 


f ude = f ude. 


C T 


In order to proof this relation, we take into account that [ and C can in a first step 
be connected by two auxiliary lines in such a way that two simply connected domains 
G, and G» with the contours C and C3, respectively, will occur. Because of 
curl w=O in G, as wellas in Go, we see in a second step that 


f udx = 0 as well as f udx = 0 
C1 C2 
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hold. This leads to 


f ude + $ ude = 0 


Ci C2 


(cf. Figure 6). 


L» 








Figure 6: Auxiliary step in the computation of the buoyancy generated by plane 
potential flows 


We realize in a third step that the integrations back and forth along each of the auxil- 
iary lines extinguish each other such that 


f ude + $ ude = f ude + $ ude 
C 


1 C2 a 


results. Here, —I’ denotes the contour of the rigid body if run through in negative 
direction, and this ends the proof. 


Remark: 

It has to be mentioned that the proof given is simplified because I’ does not only pass 
through the interior of the domain occupied by the fluid but is part of its boundary. As 
a matter of fact, one really has first, roughly speaking, to investigate the case where 
T is replaced by a line Tz of distance € from IT and passing only through the fluid, 
and in a next step one has to study the limit situation £ — 0. 


Let [ now be parametrized by r=r(r), O<7<T. Then 


T T 
Z = f (abet te) a = [piirit (1.42) 
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holds, and r = ( . ) is a vector tangential to the curve I` at the point 
(x(T), y(T)) such that 


1 ; 
n= —— 4 (1.43) 
[2 + y2 =p 
is anormal unit vector of I at this point. 


According to (1.13), we assume that the velocity of the flow at the surface of the rigid 
body is tangential to this surface: 


(u,n) = 0 along T, 


Le. Uy Y = u2. From this and with (1.42), the circulation becomes 


T 


pe feniu aaa ae: 
0 
Le. 
Z = fuo dz, (1.44) 
T 
with 
w(z) = u(x, y) = ius(x, y) z (1.45) 


w(z) is a holomorphic function in the whole domain outside the airfoil because (1.39) 
yields 


Ont = Oxy = Oya -m Oy (—u2) ; 
and because (1.28) leads to 
Oyu = Oxy E Oya = —ðr(—u2) , 


i.e. the Cauchy—Riemann equations for the function w(z) are fulfilled. 


Hence, also with respect to the computation of the circulation, it is allowed to integrate 
along C instead of T : 


Z= f w(z) dz. (1.46) 
C 
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Here, we choose C in such a way that it lies in the annulus between two concentric 
circles where K, of radius r surrounds the airfoil and where Kp is a circle with 
sufficiently great arbitrary radius R > r. Without loss of generality, we assume the 
center of the circles to be the origin of the former (x, y)—plane which now became the 
complex z-plane. 


w(z) can be represented in the annulus — and therefore particularly for every z € C — 


by a Laurent series? around the center z = 0 (cf. Figure 7): 
+00 
w(z) = 5 ay z” 
v=—o0 








Figure 7: Annulus around an airfoil 


We know by experience that the fluid flow is often influenced by the rigid body only 
in a certain neighbourhood of the contour; the flow around a ship that crosses a calm 
lake with constant velocity gives an example. From this point of view we are going to 
assume that the velocity u = u(x, y) is constant for |z| — co: 


lim u(z,y) = Ga (1.47) 


|z|00 U2,00 
with constant components U1,oo and u2,» such that 


lim w(z) = Wo — iuz o. 
|z| >œ 





>2Dierre Alphonse Laurent (1813 - 1854); Le Havre 
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Because the coefficients a, (v = 0, +1, +2, ...) are constant numbers which do 
not depend on R, and because R is allowed to tend to infinity, we find for |z| — oo 
that all the a, with positive index have to vanish. 





Hence, the result 


lim w(z) = ao = Uo — iUz. (1.48) 
|z|+00 


follows. 


On the other hand, Cauchy’s formula 


1 
eS ay fous, 
C 


leads to 


Z = 2min- (1.49) 


The (twodimensional) force K to be calculated caused by the flow and acting on the 
rigid body turns out to become 


K = ee 7 [rma (ix = Eee) l (1.50) 


where n is the unit normal vector from (1.43), p denotes the pressure along the surface 
of the body caused by the fluid, s measures the arc length along I beginning at an 
arbitrary point of it, and where L is the total length of I’. Let us study this force 
separately from the other exterior forces acting on the body, i.e. let us assume that the 
sum of these other forces vanishes. This leads to a constant potential U in Bernoulli’s 
equation (1.31) and therefore to a constant total pressure 





p 
po := 5 lul? + p. 


42 


L 
Because I’ is a closed contour, f nds = 0”. Therefore, 
0 


L 

K = — | {po -£ lul?} nas = $ f rnas 
0 

ky = 


Va? +y? dr 


The temporary introduction of the complex number 


T 
1 d 
K = -$ f —— |lu? s É ar. 
0 


k = K+ikK 


yields 


Ma eng 


The complex number within the paranthesises on the righthand side of 
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(1.51) 


lul? (@-a9) Zar = -£ sj lul? (2-49) dr 


lul’ (t -iġ) = (up +u) (t-iġ) = (uı— iuz) (u + tue) (t —iġ) 





= w(z) (u1 t+ uz2ğġ+iuzģt — iu ğ) 


is not really complex because the imaginary part vanishes (cf. (1.13)). It can therefore 


be replaced by its conjugate complex number: 


lull? (t +iġ) = w(2) (ut +uzġ-—iuzi+iu y) 


= w(z) (u; — iuz) (t +iġ) = w?(z)ż. 





— = yt? + y? , hence 


fes = | see (2) Se = (205720) 


1.3 Potential Flows and (Dynamic) Buoyancy 43 


This leads to 


2 
2 


T 
[ Po \)żdr = a i w(z) dz. 
0 T 


Because w is holomorph, also w? is a holomorphic function. Hence, 
2 p 2 
k = E i w(z) dz. (1.52) 
C 


But within the annulus we have 





a a 
w%(2) = (a + +H ) (ag + S$ + +. .) 
ag aâ—1 A 2 A 3 
= 2 = Se es 
ee ra eer aa 
with certain coefficients A, (v = —2,—3,...). 


Cauchy’s residuum formula and (1.52) therefore lead to 
ps ck 2a9 a1 + ri. 
2 
Thus, (1.48) and (1.49) yield 


k = —p (ui œ —tU2,00) Z, 


and because the circulation Z is real (cf. (1.27)), the comparison of the real and 
imaginary parts of (1.51) results in 


Kı = PU2,00 Z 
(1.53) 
Ko = —p U1,00 Z. 


If particularly u20 = 0 but u1,œ #0, i.e. if the undisturbed flow is parallel to the 
x-axis, and if Z #0, we find a lift Kə Æ 0, i.e. there is a force acting on the rigid 
body perpendicular to the direction of the flow. 


By suitable construction of wings, the airfoils can lead to Z < 0 so that an air- 
craft can start against its weight, a hydrofoil can lift in the water etc. Of course, 
also Archimedes’ static buoyancy, given by (1.62), has additionally to be taken into 
account. 
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Remark: 
The formulas (1.53) are called Kutta-Zhukovsky buoyancy formulas.”* 


Remark: 

Whereas the first equation in (1.53) describes the buoyancy at least qualitatively in a 
correct manner, the result Kı = 0 contradicts real experiences: Also in case of an 
incompressible irrotational stationary flow —as more or less realized by calmly running 
little rivers- the flow leads to some force affecting the rigid body —e.g. a bridge pier- 
parrallel to the flow. Of course, this contradiction results from one of our idealisations, 
namely from the assumption of an inviscid and therefore frictionless fluid. 

In order to understand what really happens parallel to the flow, we have to reduce these 
idealizations. This will later be done by the transition from the Euler equations to the 
so-called Navier-Stokes equations and from (1.13) to the so-called no-slip condition 


u = 0 along I. (1.54) 


(1.54) expresses the idea that moving viscous fluids leave a monomolecular non- 
moving layer along impermeable walls of solid bodies because of adhesion. Thus, 
friction along such surfaces does not mean friction between the fluid and the solid 
material but always just friction between fluid particles. 

Of course, there are also other boundary conditions where partial slip on the surfaces 
occur, e.g. in case of rarefied gas flow, porous walls etc. The tangential component of 
the flow is then proportional to the local shear stress. 


If (1.33) is reduced to the case of an irrotational, stationary, plane flow, i.e. u3 = 
0, Oyu, = ruz, if there are no exterior forces, and if (u1 , u2) is for convenience 
replaced by (u , v), we get 


(1 E (=)’) Oyu + (1 5 (°) Oyv — Sr . (ðu + əðzv) = 0 . (1.55) 


The u-v-plane is called the hodograph plane. 


We assume the equations 


u = u(x, y) 
v = v(x, y) 
to be invertible such that x and y can be expressed by u and v, i.e. 
_ | ôru Oyu 
Di sv Oyv foe 








Martin Wilhelm Kutta (1867 - 1927); Stuttgart Nikolai Jegorowitsch Zhukovsky (1847 - 1921); 
Moscow 
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In other words, it will be assumed that the vectors 0,u and Oyu are linearly indepen- 
dent. 


This leads immediately to 
ru = D vy , Oyu = —D y£ , Orv = —D Ouy , Oyv = D Oye . 


The nonlinear equation (1.55) together with the equation 0,u = O,v for the functions 
u(x, y) and v(x, y) of the irrotational flow can so -after division by D- be transformed 
into the linear equations 


v? u? uv 
Oyx — uy = 0 (1.56) 
for the functions z(u, v) and y(u, v). 


The transition from the original equations to the linear equations (1.56) is called 
method of hodographs and corresponds to the Legendre transformation in the the- 
ory of partial differential equations. 


In simply connected domains, x(u) can be derived from a potential O, ie. £ = 
Ou, y = Ov, because of 


curl æ(u) = (0, 0, Quy — Ax)’ = 0. 
(1.56) therefore yields 


2 
(1 = =) Puu + (1- =) Aw + 257 OwO = 0. (1.57) 
Cc C2 


If polar coordinates (w = ||w|| , œ) are used in the hodograph plane, i.e. 


u=wcosa , v=wsina , 
1, 

Oy = COSa Oy — — sina Oy , 
w 


1 
Oy = sina Ow + — cosa ôa , 
w 


(1.57) turns to the so-called hodograph equation 


On +3 (1-5 Sr) Ano +2 (1- g) ©=0 (1.58) 
w 


which does not explicitely contain a. 


46 1 Ideal Fluids 


If we try to solve (1.58) by the ansatz 
O(w,a) = g(w)sin(ma) or O(w,a) = g(w)cos(ma) (m€ R), (1.59) 


we receive for the unknown function g(w) the ordinary differential equation 
1 w? m? w? 
g'(w)+— (1 — =) g'(w)- — (1 — =) g(w) =0. (1.60) 


Solutions of the hodograph equation of type (1.59) are called Chapligin solutions”. 


Let gm(w) be a solution of (1.60) that belongs to a particular m and let Om(w, a) 
be the solution of (1.59) that corresponds to this solution. 


Examples: 
m = 0 leads to 


go” 1 w 
= pt at aed 
go w ĉ 
hence 
j c0 lu? 
go =——e? 2 . 
w 


pi w2 
The power series of e? & converges uniformly for all values of w. Integration can 
therefore be done term by term and yields 


ee) wy 2v 
go(w) = c finu + 5 st + c0 . 
v=1 


Here, cl) and cl) are arbitrary constants. 


Analogously form = 1: 


_ (1) aj 1 wlaw Ve wet) 
ae (art geet en 


Because (1.58) is linear and homogeneous, also all linear combinations of particular 
Chapligin solutions solve the equation (1.58). The coefficients of the expansion as 
well as the constants lm) and cv") have to be chosen in such a way that the expansion 


fits the given situation at least approximately. 





°5C.A. Chapligin: Sci. Ann. Univ. Moscow. Math. Phys. 21 (1904), 1-121 
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1.4 Motionless Fluids and Sound Propagation 


Obviously, because viscosity does only play a role in moving fluids, the results of this 
section are valid for real fluids, too. 


Let us —in a first step— consider the case of constant density p, approximately realized 
in liquids. 


In this situation, (1.26) leads for motionless fluids, i.e. for u = OQ, to the so-called 
hydrostatic equation 


and (1.31) becomes 
U+ L= const . 
p 


Let us assume that we do not yet know how the free surface of a motionless liquid 
behaves if only the force of gravity and a constant exterior (e.g. atmospheric) pressure 
po affects this liquid. Hence, the force per mass unit is given by 


k= (0, 0, gi (g = acceleration of gravity) , 


such that U = g z + const, i.e. pg z + p= const. Particularly, if (x, y, zo) isa 
point of the surface, pgz + p = pgzo + po or 


pg(z— zo) = —(p — po) = -Ô (1.61) 


holds. Here, p is the overpressure inside the liquid compared with the exterior pres- 
sure. 
Because of 
const — po 
20 = ——; 
Pg 


zo is constant, thus independent of (x, y). In other words, the surface of the liquid is 
plane or —more precisely— parallel to the earth surface. If h = zp — z is the height 
of an arbitrarily shaped liquid column, (1.61) gives 


pghq=pq, 


where q is the base area of the column. 
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The force affecting this base is given by the righthand side of the equation and does 
not depend on the form of the column whereas the lefthand side gives the weight of a 
cylindrical column of the fluid of the same base area and of the same height. 


This phenomenon is called the hydrostatic paradox. 











Figure 8: Hydrostatic paradox 


We are now going to dip a solid body of volume V and of surface F' into a non-moving 
liquid. 


Obviously, an overall force K with 
K=(k,, Kə, Ky)" =- [pndo+G 
F 


affects the body where p is the interior pressure of the liquid. G is the weight of 
the body, and denotes the outward directed normal unit vector at the points of the 
surface. Because of G = (0, 0, —G)", (1.61) yields 


Ki = pg f(e- x- Zm do = pg | (a,n) do 
F 


F 


T 
with a := (z — zo — zsh , 0, 0) . The divergence theorem therefore leads to 
Pg 


Kı=pg [ävaav=pg [E2 wv = 0. 
V V 


Analogously: Kə = 0. K3 turns out to be given by 


Ks=pg |E ave = pgV -G, 
V 
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such that 
K = (0, 0, pgV—G)!;. (1.62) 


p is the density of the liquid (!) so that pg V is the weight of the particular part of 
the liquid which is displaced by the solid body. Hence, the body is affected by a 
force against its weight direction, and this (static) buoyancy equals the weight of the 
displaced quantity of the liquid (Archimedes’ principle)”°. 


In order to study the sound propagation in a fluid, we assume that the fluid does not 
move from a macroscopic point of view. Moreover, we will only consider such sound 
effects which are due to very small disturbances of the density and of the pressure 
within the fluid. So, we do no longer assume the density to be constant but to be of 
only small variations. 


Let ô and p be the averages of the density and of the pressure, respectively. Only these 
averages are expected to be constant quantities. 


The disturbances of the density will be expressed by 
p= p (1+ o(@,t)) 
with |o| <1 and with only small spatial derivatives of o”. 
The continuity equation (1.4) then becomes 
port+pdiv (l+o)u) =0, 
thus, after division by # and taking the assumptions on ø into account, 
o.tdivu=0. (1.64) 
Experiments show that the propagation of sound waves occurs more or less in an 


adiabatic way, i.e. without gain or loss of heat. The equation of state to be taken into 
account therefore reads as 


pp = const? (1.65) 


such that 


SORSA 


2 Achimedes (about 220 BC); Syracuse 

°7Our results do not necessarily hold if the propagation of great variations of the pressure or of the 
density are under consideration as they can occur in the case of detonations. 

8 from (1.9) 
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Because ø is small, 
(lto)’el1l+yo. 
Hence 


p=pl+yo) , 
so that 

1 BY 

= Vp = ———. Vo. (1.66) 

p (1+0) p 
In our model of small disturbances, the velocity of the fluid particles and its spatial 
derivatives are so small that higher order terms of these quantities can be neglected 
compared with first order terms. In other words, the convection term in (1.6) can be 
neglected. Also exterior forces do not play a role in our context. 


If then (1.66) is put into (1.6), we find 


nee EN 
pP 


where 1 + o was approximated by 1. 


Forming the divergence of this term followed by a change of the sequence of the time 
derivative and the spatial derivatives, we end up with 


This result can be compared with (1.64) when differentiated with respect to t. This 
comparisons yields 


= 2 hes 
P 


a wave equation for ø. It shows that the sound waves propagate within the fluid with 
the velocity 


> 


AEN (1.67) 
p 


Definition: ¢ is called the local sound speed. 


Remark: Because of (1.65), ĉ can also be represented by 


ê = Ai (1.68) 


2 Weak Solutions of Conservation Laws 


2.1 Generalization of what will be called a Solution 


As already announced, we are going to discuss in this chapter the introduction of a 
suitable globalization of the definition of a conservation law solution. And we do 
already know that such a solution can’t be expected to be a smooth function because 
we had already realized that discontinuities can arise even from smooth initial states. 
Moreover, discontinuities which arose at an instant tg can move with respect to space 
and time as will be seen later. Of course, this also holds if the initial state already 
shows one or more discontinuities. 


On the other hand, discontinuities which existed at an instant tọ can disappear for 
t > to as far as a suitable definition of global solutions will be used. Examples will be 
presented later. All these phenomena can occur in scalar problems as well as in case 
of systems of conservation laws. 


We restrict ourselves within this chapter to only one space variable but allow the 
solution V = (vy,-:- Um)” seeked for to consist of more than one component. 
Hence, we are going to treat the problem 


OV+0:f(V)=0 on Q={(x,t)|rEeR,t>0}, 
(2.1) 
V(z,0) = Vo(2) , 


where we will assume for the time being that f is at least one time differentiable. 


Example: 
Another example of this kind —besides (1.14)— is the system of the so-called shallow 
water equations. 


These equations describe the shape and the velocity of surface waves as far as shallow 
liquids are concerned. Herewith, we consider an idealized situation, namely waves 
of an inviscid incompressible liquid travelling only along a canal of low depth whose 
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direction coincides with the x-axis whereas the flow of liquid particles into z-direction 
can be neglected. This idealization seems to be acceptable if the amplitudes of the 
surface waves are small compared with the wave length. The width b of the canal is 
assumed to be constant and we ask for the hight h(x,t) of the liquid surface over the 
bottom of the canal (cf. Figure 9). 

















x 


Figure 9: Shallow water flow along a canal 


The mass of the fluid within a segment x < € < x + A z of the canal at the instant t 
is given by 
atAaz 
m(t)=bp | MEDE. 
T 


With (1.3), if y is particularly chosen as p = b ph, conservation of mass leads to 


d d atAaz +AT 
Em) = ber | Mende=bo | {dh aluh) ag 
=0 V(x, Ax), 


This replaces the continuity equation in this particular situation. 


The momentum of the mass m(t) reads 


+AT 


bp / ulé,t) h(E, t) dé . 
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The forces per length unit occuring from the hydrostatic pressure, namely 


h 
h2 
bp = fo —z)pgdz=pg Ps (g : acceleration of gravity) , 
0 
act on the liquid in the segment under consideration. The principle of conservation of 
momentum then leads by means of (1.3), now with p = bpuh, to 


A(phu) + On(phu*) + Ap = 0, 


Le. to 


h2 
lhu) + Op(hu® + 5a) = 0% 


This equation turns out to be Euler’s equation for our example. 


If the continuity equation as well as the Euler equation is multiplied by g, if O,(hu) 
is replaced by ðh u + h ðu, if a variable y(x, t) := g h(x, t) is introduced, if the 
continuity equation is again taken into account, and finally, if the second equation is 
divided by w, we find the system of shallow water equations 


Y uw ` 
TEES 


Obviously, this is in fact a particular system of the type (1.2). 


In order to establish now a suitable generalized definition of a weak solution, let us 
temporarily assume that there is a smooth solution of problem (2.1) on Q. Put this 
solution into (2.1), then multiply (2.1) by an arbitrary test function ® € C(Q). 
Here, C (Q) is defined to be the set of all functions continuously differentiable on Q 
and with compact support. 


Particularly, each of the functions ® vanishes on the boundary of its support, possibly 
with the exception of such parts of the boundary which belong to the x axis. 


Now integrate by parts over Q. As a matter of fact, this is only an integration over the 
compact support of ®, hence over a closed bounded region. 


This integration leads to”? 
+00 
fiva P+Ff(V) 3P} dæ, )+ f Volx) ®(2,0)dx =0 V@OE CHO), (2.2) 
Q —oo 


>°Tt should be noted that this reformulation of the problem uses the conservation form of the differen- 
tial equation and can’t be generally applied to arbitrary quasilinear partial differential equations. 
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Obviously, the left hand side of (2.2) does no longer make sense only for sufficiently 
smooth functions V but —for instance— also for all V € L'°°(Q), i.e. for all vector 
valued functions V, whose components can be integrated over compact subsets of Q 
in the sense of Lebesgue. 


Let us now forget the way we got (2.2), and let us —vice versa— ask for functions 
V € L (Q) which fulfill (2.2). These functions will be called weak solutions of the 
original problem (2.1), and these solutions are no longer necessarily smooth. 


On the other hand, our motivation of (2.2) shows that every smooth solution is also a 
weak solution. 


Vice versa, every weak solution V which is a smooth function in a neighbourhood 
D(P)) of a point Po = (xo, to), does there not only fulfill (2.2) but also (2.1) in 
the classic understanding. Indeed, let © € Cj(Q) be an arbitrary test function with 
support D(Po). Then 





f vaT + IV ordet + f Vo2e@.0) Heo 
D(Po) -o 


holds and leads to 


{OV -D — ô (VD) + dz F(V): $- Oz (F(V)- ®)} d(x,t) — 
D(Po) 


By the divergence theorem, the last equation can be written as 


f (wta f(V)} Dde) — I (((MY)).n)) rope) 


D(Po) OD(Po) 


(i=1,--,m), 


where n denotes the outward normal unit vector on the boundary 0D(Po) of D(Po) 
and where (-, :} means the standard R?-scalar product. 
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If OD(Po) does not contain points of the x-axis, both boundary integrals vanish 
because we have ® = 0 along OD(Pp). Otherwise, the first boundary integral turns 
into an integral over the particular part of the boundary that belongs to the x-axis 
because ® vanishes along the rest of this boundary. In each of these cases we find 


i (akvi + On fi(V)} ® dz, t) 


D(Po) 
[MEP)C)em}oonen 


i.e., for i = 1,--- ,m, 


{OAV + 0, f(V)} ®(z, t) a(x,t)+ f {V(e,0)-Volx)}®(e,0) dx = 0. 
D(Po) =00 


Since ® was arbitrary and V was assumed to be smooth, this can only hold if 
OV + 3 f(V)=0 


on D(Po), hence particularly at Po, and if V (xo,0) = Vo(xo). But because Po 
was chosen arbitrarily in its neighbourhood, our assertion follows. 


This concept of a weak solution of (2.1) leads to the advantage that the set of admis- 
sible solutions can sbe extended considerably. In particular, discontinuous solutions 
can be admitted as far as the Lebesgue integrability will not be disturbed. 


But there is also a remarkable disadvantage: The uniqueness of the solution as it was 
at least guaranteed locally —together with its existence— gets lost. Often there exists 
more than one weak solution so that a new question arises: By which procedure can 
the particular weak solution answering the originally formulated real world problem 
be picked out of the set of all the weak solutions of (2.1) ? 


As a matter of fact, the solution or weak solution of a differential equation problem 
can often only be described approximately by means of a numerical procedure. Nu- 
merical procedures normally consist in discrete models of the original problem, i.e. 
in finite dimensional problems whose solutions are expected to lie in a certain neigh- 
bourhood of the unknown solution of the original problem. Particularly, we expect 
the approximate solution to converge to this unknown solution if the dimension of 
the discrete model problem tends to infinity. If convergence of the numerical solution 
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to a weak solution of the original problem can be shown, one has to ensure that this 
weak solution coincides with the relevant solution instead with one of the other weak 
solutions.°? 


For the time being, let us in the next section study an example where the transition to 
a weak formulation does really lead to a loss of uniquenes of the solution. 


2.2 Traffic Flow Example with Loss of Uniqueness 


As far as a sufficiently dense car traffic is concerned, the density p [cars/ml] as well 
as the flux f [cars/h] of macroscopic traffic flow models are approximately regarded 
as functions defined on a continuum, e.g. fora < x < b. The spatial variable x [mile] 
denotes the positions along a oneway lane which is interpreted as a onedimensional 
set, and 0 < t < oo. 


If there are no approaches or exits and no crossings along the part of the road under 
consideration, the conservation of mass, i.e. the constant number of cars, obviously 
leads within our continuous model to the demand for the validity of the continuity 
equation 


Kp + ôf = 0 (2.3) 


with 

f = vp, (2.4) 
where v = v(x, t) [ml/h] denotes the speed of the cars at position x at instant t. 
(2.3) has to be completed by an initial condition 

p(x,0) = polz) (2.5) 


with a given function pọ as well as by a spatial boundary condition as far as the part of 
the road under consideration can’t be regarded as a lane of infinite length. Moreover, 
it should be assumed that —hopefully— friction between cars does not occur. 


As by experience, p depends above all explicitely on the speed of the cars: 


p = plu). (2.6) 


30Of course, convergence is a basic requirement and must be accompanied by other important proper- 
ties of the numerical procedure in order to make it an applicable method. 
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It is an empirically given decreasing and strictly monotone function. Hence, we can 
also consider v to be a function of p. 


Let vy be the mean value of the individual maximal car velocities if the road in front 
of a driver is more or less empty. This value is assumed to be a constant, e.g. vp = 80 
ml/h. 


Obviously, the initial density 
po(x) = plv(#,0)) 
is known as soon as the initial speed distribution 
v(x,0) = volz) (2.7) 
is given. 
Instead of p(v) or of v(p), often the graph of 
f = vlo): p = flo) (2.8) 


is given empirically and turns out to be a strictly concave function called the funda- 
mental diagram. 


Thus, this traffic flow model leads to an initial value problem for a scalar conservation 
law of type (2.1), where V consists of the only one component p such that f(V) = 


Fp) 


Let p* denote the jam concentration, i.e. the density of cars within a standing con- 
gestion. Also this value is assumed to be a constant similar to vy, e.g. p* = 400 
cars/ml. 


A most simple f(p) model, later used in order to explain exemplarily certain facts and 
connections, was given by Greenshields: 


f(p) = vg: p- (1-4). (2.9) 
Thus, the flux vanishes or —in other words— there is (almost) no traffic, if the traffic 
density vanishes or if the road users stand within a traffic jam. 

aad 


(2.9) is aconcave parabola with its maximum at ( eae 


) and leads together with 
(2.4) to the relation 


D= p (1-2), 
Uf 
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Hence, in case of Greenshields’ model, the graph of p(v) turns out to be a straight 
line. 


The more general model for arbitrary fluxes is due to Lighthill & Whitham?!. 
By means of Greenshields’ model we are going to start from the top of a traffic jam, 


i.e. 


(2.10) 
0 for r> 0 


for to = 0. The top of the traffic jam (red traffic light) is located at x = 0, and we 
identify the forward direction of the traffic with the positive direction of the x-axis. 


Remark: 


Conservation law problems whose initial values are constant on the left hand side of a 
discontinuity and constant also on its right hand side are called Riemann-problems””. 


With (2.9) and (1.17), the characteristics are given by 


1 —t for xo <0 
t for xo = 0 











xT 


Figure 10: Characteristics in case of the dissolution of a traffic jam 


31M_J. Lighthill and C.B. Whitham: Proc. Roy. Soc. A 229 (1955), 317-345 
Bernhard Riemann (1826 - 1866); Göttingen 
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Obviously, the characteristics do not intersect for t > 0 so that we find immediately 
p* for xz < —vpt 
p(z,t) = 
0 for x> vet . 


Thus, the question is open how the solution will look like within the hatched region 
of Figure 10, i.e. for 


SUF < x < vst, t>0. 


It can easily be verified that 


* 


pa p fr «<0 
piso) = { 0 foi? ai Ve 20 (2.11) 


is a weak solution of 
R _ 
Op + fy + Ox (oa - 4) = 0, 


i.e. of the Greenshields type of (2.3) together with the initial condition (2.10), because 
the left hand side of (2.2) if realized for our particular flux f leads to 


i J [2:81 tva- D] dlet) + f O(c, 0) pol) ae 


p* Ms 
œ 0 0 
=f f pw aSd(a,t) + f (2,0) ~ dz 
0 =œ =00 


0 oo 0 
pe dag dt +(e, 0)} dx = p* f {®(x,0) — ®(x,0)} dx =0. 


—oo 


In other words, (2.1) is a solution which keeps its initial state: The drivers do not start 
though the traffic light switches from red to green. This solution describes a special 
situation where the state along the left hand boundary of the hatched region of Figure 
10, i.e. along the characteristic t = ——, is connected with the right hand boundary 
v 

x : : ae ‘i 
t = — across a line of discontinuity. 

UF 
Of course, drivers do not behave this way, so that (2.11) is certainly an unreasonable, 
hence irrelevant solution. As a matter of fact, the drivers ride impulse was not yet 
introduced into the model. 
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But there is also another weak solution of our problem, namely 


* 


p for x< HUE, t>0 


„Uft— r 


polz, t) := rae 
uy 


for eta oe <vet, t>0 , (2.12) 


0 for z>ust, t20 


as can also easily be veryfied by elementary integration procedures. Within this solu- 
tion the states along the left- and the righthand sides of the hatched region of Figure 
10 are continuously connected. Herewith, the density decreases so that this type of a 
solution is called a rarefaction wave. 


The head of the traffic jam, i.e the the front end of the still standing chain of cars, 
withdraws against the direction of the movement of the traffic with speed —v+¢ whereas 
the front cars move with speed +vf. Obviously, (2.12) coincides with our experiences 
within the bounds of Greenshields’ idealizations. It therefore has to be regarded as the 
particular weak solution that fits the circumstances of the given real situation.’ 


Of course, the loss of uniqueness can also occur if systems of conservation laws in- 
stead of scalar situations are concerned. 


As a further application of Greenshields’ traffic flow model, let us consider the prop- 
agation of the end of a traffic jam opposite to the direction of the movement. 


We start from an initial situation at the instant t = 0 where the traffic shows maximal 
flux on the left of x = 0, ie. p = as whereas the cars on the right of = 1 


did already stop, i.e. = p*, and where we have a continuous linear transition in 
between: 


Z for z<0 
polz) = Sa +2) for O0<a@< 1 ao) 
p* for g> 


33This is also an example of a situation mentioned earlier, namely that a solution being continuous for 
t > 0 can develop from a discontinuous initial function. 
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The characteristics then turn out to be 





£T = T0 for zo <0 
1 x 
t=— (1 = =) for 0<% <1 l (2.14) 
Uf XO 
pene for to >1 
Uf 
t 
1 
vf 
x 


Figure 11: Characteristics at the creation of a traffic jam 


1 
This solution p(x, t) keeps its continuity for the time 0 < t < —. Then suddenly 
Te 


; Puy 1 prai ices 1 
a discontinuity occurs at t = —, and this discontinuity moves for t > — along 
v v 


a not yet known curve I. But E knows that this curve takes a course ine left 
upper half plane because the characteristics arising from the x-axis on the right of 
xo = 1 certainly intersect with the charactics arising from rg < 0. As soon as [ will 
be known, the values of p given along the characteristics make possible to compute 
explicitly a solution p(x, t). It is discontinuous along I but still an element of L/°, 
thus really a weak solution. 


The drivers do not expect such a discontinuity to crop up, and this often is one of the 
reasons of a rear-end collision. 


How T can be determined will be shown by one of the examples of the next section. 


62 2 Weak Solutions of Conservation Laws 


2.3 The Rankine-Hugoniot Condition 


We are going to treat the particular case where a certain discontinuity of a weak solu- 
tion V of (2.1) moves continuously along the x-axis. In other words, the set of points 
of Q where these discontinuities occur is a continuous curve I in Q. Here, it can hap- 
pen that all the components of V along this way are discontinuous or that this only 
holds for some of them. In gas dynamics, often the density is particularly concerned 
and forms a shock. 


Let P) = (£o, to) € I with to > 0 and assume a bounded neighbourhood D( Po) C Q 
of P to exist which does not contain points of the x-axis and where V is smooth 
outside I+, 


tah 


to] 











= 
Xo X 


Figure 12: Sketch to explain the proof of the Rankine-Hugoniot condition 


We decompose D(Po) into the two parts Dı and Do (cf. Figure 12) each of which we 
are going to treat as closed domains. 


Taking the smoothness of V in the open domains Dı and Dg into account, we find 
on D; and also on D3 for a particular ® € C4 (Q) whose compact support coincides 


34 The considerations to be formulated here can easily be generalized in order to fit also more com- 
plicated situations. 
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with D(Po) = Dı U Dz the relation 


z [ivaa aa 1 = J Voix) 2,0) de 








Q —0oo 
ae I a(V ®) + 3l f(V) ®)] d(a,t) + f [A.V © + 3 f(V) è] d(z, t) 
Dı Dı 
z / &(V 2) + 3l f(V) ®)] d(x,t) + i [V1 ® + def(V) 8] d(a,t) 
Də D2 


Together with the validity of (2.1) in Dı and in Do, 


law ®)+0x(f(V) ®)] d(x, t) + fow ®)+0r(f(V) ®)] d(x, t) = 0. 


Dı D2 


The divergence theorem applied to each of the components of each of both integrals 
then leads to 


J iVa) dæ- (HV) a) at} + J iVa) dæ- (FV) Zik 
OD, OD2 


Because ® vanishes on the boundary of its support, i.e. on the boundary of D, we find 
Q2 Q2 
[iVa -EVD - fV) de- FVS) dt} = 0, 
Qı Qı 


where V p denotes the lefthand limits of V along I if for every particular fixed value 
of t the values of x approach T from the left, i.e. Ve(x, t) = V (x—0, t) for (x,t) Er. 
Analogously, V „ means the values of V if we approach T from the right. The integrals 
have to be understood as line integrals along r. 


Let 
[V] := Ve- Vr , [F] := fe- fr = fVo) - F(Vr) 


be the jumps. 
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Then, obviously, 


Q2 d 
x 
ViI— - dt = 
[tig - i} | 
1 
dx F aoe 
where T =: 0 denotes the velocity by which the discontinuity moves along the x- 


axis. Since the test function ® with support D( Po) was arbitrary and because Pp was 
arbitrarily chosen point on I’, the so-called Rankine-Hugoniot condition 


[V]6 = [f] (2.15) 


turns out to hold along [ componentwise under the conditions we had assumed to be 
fulfilled. 


(2.15) is also called jump condition, and for every weak solution V there is a partic- 
ular velocity: 6 = ê( V). 


Example: 
If (2.15) is applid to our one dimensional gas flow where 


q 
P 1 > 
vae aT Ea +p |., (2.16) 
E E+p 
p 


and ifq = pu —with the flow velocity u— is taken into account, we find from the first 
component of the jump condition 


(pe — pr) Ô = qe — qr = pe ug — Pr Ur 


pe (Ô — Ue) = pr (Ô — ur) . (2.17) 
The second component leads to 


2 
A q 
(e-g)o=+-% 
Pe Pr 


+ pe — Pr, 


Pe u(Ô = ug) = Pr Ur (ô = ur) + pe — Pr, (2.18) 
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and the third component gives 





(Ee — Er) ô = pA T ye napa = ee 
Pe Pr Pe Pr 
i.e. 
Eg (ô x ug) = k, (ô ur) H pe ug — Pr ur. (2.19) 
If 
ug =U =U a (2.20) 


(2.17) is fulfilled, and (2.18) then gives pẹ = p, such that also (2.19) is realized. 


Hence, in this case, the discontinuity ’swims’ with the flow; this means that there are 
no gas particles which cross the point of discontinuity though pẹ = pr or Ey = E, do 
not necessarily hold. 


The situation described by 


Pr =U Pe» Pe # Pr (2.21) 
is called contact discontinuity. 


Otherwise, ue A ô, uy # ô show that there are gas particles which cross the discon- 
tinuity because (2.17) causes in this case that ô — ue and 6 — ur are both positive or 
both negative*>. This type of discontinuity is called a shock and also the curve T is 
called a shock or a shock curve. The particular side of the shock where the particles 
having already crossed the shock are situated is called the back side of the shock, the 
other part is called the front side. If, for instance, ô > ur > 0, the righthand side is 
the front side, the lefthand side is the back side. 


Assume a weak solution to be smooth outside the discontinuity curve I, hence to be 
outside [I also a solution of (2.1). The weak formulation (2.2) of the problem under 
consideration is then equivalent to the demand that the Rankine-Hugoniot condition 
is fulfilled along I’ and that simultaneously the original version of the initial value 
problem (2.1) is fulfilled, too. This equivalence follows from the already proven fact 
that weak solutions are at the same time classic solutions on domains where these 
weak solutions are smooth. 


In case of a linear system 





>We omit the unnatural case p < 0. 
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the Rankine-Hugoniot condition leads to 
[V] -ô = |[AV~,—-AV,] = A[V]. 


Thus, 0 is an eigenvalue of A as far as a discontinuity of the initial function occurs, 
i.e. as far as [V] A O turns out to hold for t = 0 and, hence, also for increasing values 
of t. This leads to 


t = Xr, foraparticular k 


since (t) = @ is the differential equation for I. 


But because of (1.18), this is also the differential equation for the k-th set of character- 
istic curves so that discontinuities move along characteristics as far as linear problems 
are concerned. 


Let us wind up this section by the determination of I occuring in Greenshields’ traffic 
flow model in case of a traffic jam increasing opposite to the flow direction. 


Because of the relation 


presented in the Greenshields model, and taking (2.15) into account, I’ has to be de- 
termined by the differential equation 


_ lob _ ( aie) 
r a a N |] E 
Pe — Pr P 


From Figure 11, the initial condition reads as 
r(+) = 0. (2.22) 


Moreover, Figure 11 shows 


1 
p=p for t> — 
UF 
along the characteristics coming from the righthand side, and this does not depend on 
the location of I’. Along the vertical characteristics coming from the left side, we find 
analogously 
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Hence pz = and therefore 


ae 
m 


Together with (2.22), the shock turns out to be the straight line 


t= 


1 
t = —(1-2z), «<0. 
Uf 








>= 
X 


Figure 13: Shock when a traffic jam increases 


Knowing I and the values of p along the characteristic curves, a weak solution of our 
traffic flow example can now immediately be constructed. Particularly, the propaga- 
tion of the end of the congestion against the direction of the cars which have not yet 
come to a halt can so be evaluated: 





p* 1 

pi for 0 < t< z(= 2x), —-o<2<0 
p(x, t) = pe a < lex < 

po eae for 0<t< z> OSt<l 

p“ otherwise . 


(2.23) 


This solution shows under the assumptions of the Greenshields idelizations and for the 
particular initial situation under consideration that the end of the chain of non-moving 
bs Ix . ; 
cars originally located at x = 1 will be located at x = 0 after — time units. From this 
Ue 
very moment on the end of the traffic jam moves with the congestion velocity $v f 
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against the direction of the cars still moving, and the transition from the density $v f 
of the moving traffic to the density p* of the set of cars which already came to a halt 
occurs very suddenly, at least from the driver’s point of view. 


With respect to the first part of our traffic flow model as well as generally, the ques- 
tion arises how in case of non-unique weak solutions of (2.1) the particular solution 
relevant for answering the real world problem that underlies the mathematical model 
can be characterized. 


We are going to treat this question in the next chapter. 


3 Entropy Conditions 


3.1 Entropy in Case of Ideal Fluids 


The loss of uniqueness caused by the installation of the concept of weak solutions 
leads to the necessity to formulate a criterion by which the physically relevant weak 
solution can be picked out of the set of all the weak solutions. We therefore ask for an 
additional constraint that characterizes this particular weak solution. 


In order to get an idea, let us look again to the important application of conservation 
laws to the description of the dynamics of ideal gases. Besides the state equation, we 
took into account the conservation of mass, the conservation of momentum and the 
conservation of energy. The conservation of energy was formulated by the first law of 
thermodynamics, but what we did not yet respect is the second one and its assertion 
on the behaviour of the state variable 


S = entropy/mass . 


Without going into details, it should only be mentioned that the entropy describes a 
measure for the probability of the existence of a particular physical state (Boltzmann 
statistics in thermodynamics”), 


37 


The Second Law of Thermodynamics”’ can be stated as follows: 


If there is a closed physical system without supply of energy from outside, every phys- 
ical process inside the system takes a course such that the entropy does not decrease*®, 


and it increases if possible*?. 


Figure 14 shows an example for the validity of the second law, namely a shock tube 
divided by a membrane into two parts one of which is filled by a gas whereas the other 


Ludwig Boltzmann (1844 - 1906); Graz, Vienna, Leipzig, Vienna 
37R.J.E. Clausius (1822 - 1888); Zurich, Würzburg, Bonn 

38e.9.: cf. R. Ansorge und Th. Sonar: ZAMM 77 (1997), 803-821 
irreversible processes 


Mathematical Models of Fluiddynamics. Rainer Ansorge 
Copyright © 2003 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim 
ISBN: 3-527-40397-3 
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membrane 


vacuum 


ORR) . 
LC arr be ee 





Figure 14: Example for increasing entropy 


one is evacuated. When the membrane will be taken away, parts of the gas will begin 
to move into the evacuated part of the tube though all the conservation principles will 
also be respected if the gas particles would keep their positions. And as a matter of 
fact, the Euler equations will indeed also be fulfilled if nothing would happen. But 
this, obviously, does not describe the physically relevant solution as everybody knows 
from experience®®. The reason is that the equidistribution of the particles over the 
whole tube is a state more probable than the initial situation. 


In other words, we expect the second law to be a model for the formulation of a 
constraint that allows to pick out of the set of weak solutions of a conservation law 
system the relevant solution, even if problems of other fields of applications such as 
problems of economy models are concerned. In any case, such a constraint will be 
called an entropy condition and the solution picked out will be called an entropy 
solution, hopefully the only one by the proof of a uniqueness theorem. 


Generalizations of the second law in order to find a constraint that does not only work 
for physical tasks are due to several authors, particularly to Oleinik and Lax who 
started treating this question in a first step for scalar conservation laws. 


To understand their ideas, let us look for the situation in gas dynamics, for convenience 
in case of a one dimensional flow. 


It was already stated in section 1.1 that each of the state variables can uniquely be 
expressed in terms of two of the other state variables by means of a so-called equation 
of state; a particular example was demonstrated by (1.9). 


One result of the theory of thermodynamics leads to the relation 


E de+pdW 


ds T 


(3.1) 


“By the way, this is an analogue to the situation where a traffic light switches from red to green but 
the drivers do not start. 
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x ~ 1 
where W denotes the specific volume W := — (volume of mass 1). 
p 


If we ask particularly for the connection between the three state variables €, S' and 
p= T the equation of state belonging to these variables leads to 


de = Ose dS + Oye dW , 
hence with (3.1) to 

TdS = dsedS + ðe dW + pdW = Ose dS + (Ope + p) dW , 
such that 

Ose = T, Owe = -p (3.2) 
follows from the linear independence of dS and dW. 
Because of 

He = Ose OS + Owe OW = TAS — poy » 


the relation 


o, 
TOS = Ope + pow = Ove = poy 
follows. Analogously: (3.3) 
Ox 
T ôS = rE — p l 


Denoting the entropy per volume by s := pS, the quantity ôs + O,(us) can be 
expressed as 


Os + slus) = Alps) + Oz((pu) - S) 


rp S + pS + ôrlpu) S + purs 


The first term on the right side vanishes because of the continuity equation, and taking 
(3.3) into account the last equation leads to 


Os + Oz(us) = £{ (a-p) + u (a-p) : (3.4) 
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The conservation law (1.7) concerning the specific total energy Æ and formulated in 
the particular situation of a one dimensional flow reads as 


or — with E = pe + p u?, reformulated by means of primitive variables — 
1 2 
Osp € + poe + 3 ıp u + pu dru 


1 
tu fapt dapet pret 5 OxP u? + pudeuh+ (pe Fw +p) O,u=0, 


2 


(O:p + Ox(pu)) (« + =) +p (a + uôrE + P Osu) 


+ pu (am ŽE + usu) = 0. (3.5) 


Also in this equation, the first term vanishes because of the continuity equation. The 
2 

second conservation law, namely Ojg + 0,(5 +p) = 0, also rewritten by means of 

primitive variables, leads to 


Or 
P (a+ a tuða) + u (p + Or(pu)) = 0, 


so that also the third term on the left side of (3.5) disappears, again by means of the 
continuity equation. Thus, after division by p, 


Oe + uzre + P Ant = 0, (3.6) 
p 
remains, so that the term inside the outer braces in formula (3.4) turns to 


0, On 
-Egu PM PIT _ = lapt Olou) = 0, (3.7) 
p pp pP p 


again because of the continuity equation. 


Hence, also s fulfills a conservation law: 
Os + O,(us) = 0, (3.8) 


and this automatically provided that the Euler equations are fulfilled. 
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But we have to take into account that this result required the differentiability properties 
to be ensured which we had used in order to derive (3.8). Thus, if the flow is smooth, 
section 1.1. shows that the total entropy 


J sd(x,y,z) 


W(t) 


of a volume W(t) arbitraryly taken out of the fluid at an arbitrary instant t stays 
constant with respect to time. This is known as an isentropic flow. 


Equation (3.8) can not be expected to hold if the smoothness assumptions are not ful- 
filled. But what can be expected, because of the second main theorem, is the validity 
of the inequality 


Os + Oz,(us) > 0 (3.9) 


in its weak form, 1.e. 


[eas palsp olds: f s(o.0)0(2,0)az <0 
Q R 
V@eECj(Q) with 6 > 0. (3.10) 


Formally and because of its inequality character, (3.9) was multiplied by nonnegative 
test functions, then integrated by parts, and finally this ‘derivation’ of (3.10) should 
be forgotten. 


Equation (3.10) is the entropy condition for the flow under consideration, and the 
‘derivation’ shows that it is indeed satisfied automatically provided that the original 
conservation laws are fulfilled in the classical understanding. This follows from the 
fact that (3.8) was a direct consequence of the Euler equations. 


If the weak solution V of the Euler equations is discontinuous, hence not necessarily 
unique, (3.10) is really an additional condition. 


Remark: 
We mention without proof that the physical entropy per volume —s = —s(V) of gas 
dynamics is strictly convex. 
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3.2 Generalization of the Entropy Condition 


It was Lax*! whose definition of an entropy solution of a general system of conser- 
vation laws (2.1) can in the case of one spatial variable and by imitation of the gas 
dynamics situation be formulated as follows: 


Definition: 

A weak solution V = (v1,:°: Um)” of the system (2.1), i.e. a solution of (2.2), 
is called an entropy solution of this system if there is a scalar and strictly convex 
function § = (V) as well as a scalar function F = F (V) belonging to Š so that in 
every domain where V is smooth, (2.1) leads automatically to the validity of 


a,S(V) +0, F(V) = 0, (3.11) 


and so that in the non-smooth situation theentropy condition 


l {5V 4) Blet) + ËV (e, 1) 2e, t)} det) 
Q 


+f S(Vo(x)) ®(2,0)dx >0 YS e CMQ), > 0 (3.12) 
R 


is respected. 
Remarks: 
— Š is often also called the entropy functional and F is called the entropy flux. 


— In gas dynamics we have S=-s. 


Sometimes we will more specifically talk of Lax entropy solutions because later 
also some other definitions will be given. 


— (3.12) is often abbreviated by the formulation: 
a5(V) + OxF(V) <0 holds weakly. 
After the choice of S, F is already more or less determined by (3.11). Especially in 
the scalar case m=1 (V =v, f = f) ,(3.11) leads for smooth solutions v to 
S (v) -ðw + F'(v)d,v = 0, 


ip Lax, in: Contributions to Nonlinear Functional Analysis (E. Zarantonello, ed.), Academic Press 
1971 
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and therefore with (2.1) to 

S'(v) {-Ocf(v)} + F(v) dev = 0, 
i.e. to 

-S (w) fw) + Fw) Ozu = 0. 


Because this must hold automatically independent of the choice of vp(x)* , the rela- 
tion 


F'(v) = S'(v) f'o) 
has to be fulfilled, 1.e. 
Vv 


Fa f Preda Const (3.13) 
0 


Just the constant is still arbitrary. 


If m > 1 and if the function V is smooth, (3.11) yields 
(v5, av) if (WF, dV) =y 


where 


By analogy with m = 1, putting (2.1) into this equation, we find 


— (v5, JFV) dV) + (vWF, dV) 


(0,V)" (vi = (J #(V))" Vv Š) ee 


0, ie. 


Again, because this relation must hold for any choice of Vo, S' and F have to be 
connected by 


VF — (Jf(V))' WS =0, (3.14) 


and this is a system of m differential equations for only two functions S, F to be 
determined that depends on the m independent variables v1,--- , Um. 





4# particularly also for initial functions generating smooth solutions 
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Hence,(3.14) does not necessarily have a solution in case of m > 2* so that the 
definition of an entropy condition used so far has to be modified suitably. We shall 
come back to this question in section 3.3. 


But if f(V) is a gradient, too, i.e. 
F(V) = Vv k(V) 


with a scalar function k(V), (3.14) turns by the particular choice of $(V) = 3(V,V) 
into the system of differential equations 


prs 5 dy, vik Uy = Ov, (> ante.) —~ ak (=1,---,m), 


v=1 v=1 


and this system can be satisfied by several solutions, e.g. by the particular entropy 
flux 


F(V) = (f(V), V) — k(V) + const. (3.15) 


A straight forward formulation of these considerations in case of m = 1 leads directly 
to the particular couple 


S(v) = La , F(v) = f oada + const , (3.16) 
0 


but —obviously— also an infinite number of other suitable couples exists. 


This fact generates a further question, namely for the entropy solution’s independence 
of the choice of the entropy functional S. 


In order to answer this question, we cite the following theorem whose proof is more or 
less just a word-to-word copy of the proof of the Rankine-Hugoniot condition taking 
(3.12) into account instead of (2.2): 


Theorem 3.1: 

Assume that there is a couple Š , F of functions belonging to a weak solution V and 
fulfilling condition (3.14); let the discontinuities of V form a curve T along which 
the Rankine-Hugoniot condition holds. Then it can be stated that V is a Lax entropy 
solution if and only if the inequality 


[S]é < [F] (3.17) 


“3 As far as gas dynamics is concerned, a suitable entropy flux belonging to S = —s exists though 
3<m<5. 
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holds where 6 is the velocity occuring in the Rankine-Hugoniot condition (2.15), i.e. 
the velocity by which the discontinuities move along the x-axis. 


As in (2.15), [5] := S(V~) — S(V,) and [F] := F(V,) — F(V;,) denote the 
heights of the jumps. 


In the scalar case and together with (2.15), (3.17) yields a particular form of the en- 
tropy condition equivalent to (3.12), namely 





BLA < [Ë], (3.18) 


i.e —by means of (3.13)— 


S(ve) — S(vr) 
Ue — Ur 


(f(ve) — f(vr)) - [sos da < 0. (3.19) 


This leads to the following theorem: 


Theorem 3.2: 
If the entropy functional S(v) is strictly convex“ , i.e. 


S"(v) >0, WER, 


if also the flux f is strictly convex, and if discontinuities do only occur along smooth 
curves I which do not intersect in the area under consideration, condition (3.19) does 
only hold if and only if along each of these curves I the jump relation 


Up > Up (3.20) 


is respected. In particular, this result shows that all entropy conditions formed by 
means of strictly convex functionals Š and their entropy fluxes F are equivalent be- 
cause these functionals do not play a role within formula (3.20). 


PROOF: 7 
The function f(v) defined for every fixed value of v, by 


FW) = FO) rey) — Fo) — f Sla) f(a)da for v # », 


flv) = U — Up ip 
0 for v = Up 





as it was up to now asumed to be 
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is differentiable, and its derivative reads 


a S'(v) = S(v) N S(vr) ; fro) = f) — f(vr) 


Fw) = VU — Ur U — Ur 
0 for v = vp. 


The terms within the braces are positive for v > v, and negative for v < vr because 
Sas well as f are strictly convex. This leads to 


fi(v) <0 for V É Ur 


so that v, is the only value for which f vanishes, and f is strictly monotonously 
decreasing for v Æ vp. Thus, 


fv) >0 forv<vu, , f(v)=0 forv=v, , flv) <0 forv >v. 
(3.21) 


Equation (3.19) implies f (ve) < 0, and because of (3.21) this is really equivalent to 


Up > Up. (3.22) 














Remark: 
The same arguments if applied to our traffic flow example with its strictly concave 
flux f(p) would similarly lead to the equivalence of the entropy condition to 


pes pr.” (3.23) 


and the interpretation of this result shows that car drivers endeavour to smooth a 
discontinuous situation as soon as possible, i.e. to come to an equals sign*®, or that 
they intent to achieve a nondecreasing development of the traffic density close to a 
discontinuity in being. In particular, no car driver stops when he notices that there 
is a traffic jam in some distance in front of him; he would rather drive to the end of 
the jam such that the density there will be increased. The end of the jam, i.e. the 
discontinuity, passes over him against his ride direction when he stops. 


The particular weak solution (2.11) of our traffic flow example does obviously not 
fulfill condition (3.23). In other words, (3.23) represents the driver’s ride impulse and 
completes the rough Lighthill- Whitham traffic flow model”. 

‘The strictly concave case can easily be transformed into a problem with strictly convex flux. 


“©The discontinuity at the top of a traffic jam is such a situation. 
“R. Ansorge: Transpn. Res. 24 B (1990), 133-143 
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There is also a geometric characterization of the entropy solution provided that the 
flux f is strictly convex: 


The Rankine-Hugoniot condition (2.15) together with the mean value theorem yields 


fve) = Fler) 


Ug — Ur = f'(up + O(ve — vy) , O< v= 1), 


and the strict convexity of f together with (3.20) then leads to 
f'(ur) < ê < f'(ve). (3.24) 


Thus, if and only if the weak solution is the entropy solution, the ascents of the char- 
acteristics (1.17) belonging to this solution cause the characteristics to run into the 
shock curves I for increasing t instead of leaving them (cf. Figure 15). 


t 
Tr 





Figure 15: Course of the characteristics in case of entropy solutions of scalar prob- 
lems as soon as shocks appear 


Olga Oleinik did also publish an entropy condition for weak solutions of scalar con- 
servation laws, and this already in 1957/8. It can be formulated as follows: 


Definition: 
A weak solution v of the scalar problem (2.1) is called an entropy solution in the sense 
of Oleinik —or an Oleinik solution— if there is a constant E > 0 so that 

v(x +a,t) — v(z,t) 2 E 


7 Va>0,Vt>0,VzeER. (3.25) 


Assume the discontinuities of a weak solution of a one-dimensional scalar conser- 
vation law to form a curve I’ such that the Rankine-Hugoniot condition (2.15) holds 





480, Oleinik: Usp. Math. Nauk. (N.S.) 12 (1957), 3-73 
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along this curve; if then the point (x,t) € Q (t > 0, fixed) moves to T from the left, 
the point (x +a,t) , a > 0, finally becomes a point right of T such that (3.25) states 


E 
v(x +a, t) — ve < 7 


If a now tends to zero, 
Ur — Ve < 0 
follows. 


Remark: 

Thus, an Oleinik solution is a Lax solution, too. If —hopefully— a uniqueness theo- 
rem for Lax entropy solutions will be available, the Oleinik solution is also uniquely 
determined and equals the Lax solution. 


This uniqueness problem for one-dimensional scalar conservation laws will be dis- 
cussed in the next section, and also the question concerning a suitable generalization 
of the definition of a Lax entropy solution to systems of conservation laws will then 
be treated. 


Remark: 
Another access to the idea of an entropy solution for scalar problems ensues from the 
wellknown uniqueness of the solution v,(x, t) of the parabolic problem 


Ow + Oz f(v) = €Ozzv, €>0 
u(x, 0) = v(x) ; 


if the righthand side of (2.1) is enriched analogously to the linear problem (1.24) by a 
diffusion term where the diffusion parameter € tends to 0°. 


3.3 Uniqueness of Entropy Solutions 


In this section we are going to investigate the uniqueness of the Lax entropy solution. 
Lax himself gave the outline of the following theorem: 


Theorem 3.3: 
Let ù and v be two Lax entropy solutions of the scalar form of the conservation law 
(2.1) which are piecewise continuous for every fixed t > 0, non-smooth only along 





N.N. Kuznetsov: USSR Comput. Math. Math.Phys. 16 (1976), 105-119 
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certain curves I, and belonging to initial functions vg and vo, respectively. The flux 
f is assumed to be strictly convex. Then, the distance 


0.8) = 06, Dll 


is monotonously decreasing with respect to t, where the L,-norm has to be understood 
50 


as taken over an arbitrary but fixed interval I of the x-axis°”. 

PROOF: Decompose for every fixed ¢ the interval J into intervals (x), 2,41) 
(v = 0, +1, +2,---) in such way that 0(2,t)—v(a,t) has the sign (—1)” uniformly 
on (£v, £v4+1). A decomposition of this type exists because & — v was assumed 
to be piecewise continuous with respect to x. Evidently, x, depends on t: x, = 
zult) (v=0,+1,+2,---). Thus, 


























ry+41(t) 
OC.) -= vh Ona = 5 (—1)” / [o(a, t) — v(a,t)] dz. (3.26) 
xy (t) 


v and v are piecewise smooth on (x(t), v,+1(t)), and differentiation with respect 
to t yields 


£ 6-8) = ve Dlan = do SOROR) — v(zv+1 (t), t)] 2,41 (t) 


V 


£v41(t) 
+ [Oxo — Ov (a, t)] do), 
ty (t) 
With (2.1), 
Tv+1(t) tv+1(t) 
O(a, t) d£ = — i Ox f (v(x, t)) dx 
xy (t) xy(t) 


= —flo(arsi(4),t)) + Fee), t) 
follows; analogously for v. 
Particularly, if then õo = vo in Li, i.e. ||@(-,0) — v(-,0)||;,( = 0 for every interval J, 


lõ t) — v6, tlla) = 0 follows for all t > 0, i.e. © = v in Ly at every time level ¢, and this 
means uniqueness. 
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Hence, 


= [a(ax,(t), t) = v(£y(t), t)] x(t) (3.27) 


+ [FED = Fole) 


We have to distinguish between two different situations: 


i) The function ©(-, to) — v(-, to) is continuous at x,(to) 5!. Then this function 
vanishes at x, (to), i.e. 
0(x_(to), to) = u(x,(to), to). (3.28) 


ii) One of the functions (x, to) , v(x, to) is discontinuous at £, (to). 


We therefore need to sum up on the righthand side of (3.27) only for the particular 
numbers v characterized by the fact that (x, (t), t) or (w,+41(t),t) lie on curves of 
discontinuities of v or of v. 


Let us at first discuss the case that (x ,+1(t),t) lies on a curve of discontinuities 
of ù whereas v is smooth at this position. The end point z,+1 of the interval 
(x_(t), £v+1(t)) has then to be looked upon as x)41 — 0, i.e. 0(vy+41(t),t) = ve. If, 
for example, ù — v > 0 on (zy, 2,41), hence 0 — v < 0 on (2141, @) +2), then v is 
even and 


Ug > Ve = Up > Up (3.29) 
because x14 has on (4,41, £v+2) to be looked upon as 7,41 + 0. 


This leads to 
(-1)" {flava Pen oe 
GG d OS) Siena, yi} 





5 for the arbitrary fixed value of to under consideration 
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= (ŭe — ve) B41 (t) — flu) + Fre). 


Because of the smoothness of v at (xy41(t), t), £ = x) 41(t) can only be a point of the 
shock curve of v. The Rankine-Hugoniot condition (2.15) therefore makes it possible 
to estimate the righthand side of the last equation by 


(ŭe — jp) ee) — f (ve) + f(ve) 


20D- Her) _ A-A] <o 


’ 


where it was taken into account that (3.29) yields vg — vg > 0, but that 


fe) — Fler) _ fe) — F(ve) < 


=; = 0. 
Ue — Ur Ue — Ve 


Here, the last inequality results from v, < ve < õe, also following from (3.29), as 
well as from the strict convexity of f (cf. Fig. 16). 


fi 








cT 





u, vı uj 


Figure 16: Sketch concerning Lax’ uniqueness theorem 


A similar result follows if (7,+1(t), t) is a point of a curve which is a set of disconti- 
nuities of ù as well as of v; in this situation, the inequality 


Ùp > Ve > Up > Up 


has to be taken into account instead of (3.29), and it is this part of the proof where the 
jump condition vg > vr applies to. 


The same reflections work for the values of õ and v at the position (x(t), t), namely 
disappearance or non-positivity of the terms. Here, xy (t) = x, has to be observed in 
case of a discontinuity. 
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The case 0 < v on (£y, £,+41) can be treated in the same way. So, indeed, 


Blt) — vl S 9 











follows. 





As already emphasized, in case of systems of conservation laws of more than two 
equations, pairs (5, F ) of an entropy functional Š and an entropy flux F connected 
with eachother by (3.14) and appropriate to distinguish an entropy solution from other 
weak solutions do not necessarily exist. 


We continue treating only one-dimensional problems but we try to formulate a suitable 
analog to (3.24) as far as problems are concerned for which the convexity property of 
the flux can also suitably be generalized. 


Let N C R” bea subset of m-dimensional vectors and assume the system (2.1) 
to be strictly hyperbolic for all V = (v1, v2,- Um)” € N, i.e. the eigenvalues 
Ai(V) of the Jacobian J f(V) are real and different from eachother for all V € N. 
Moreover, let f(V) € C?(N) and 


Ail V) < Aa(V) <- < Am(V). (3.30) 


The eigenvalues keep this order for all V € N since otherwise a vector V € N 
with 


Mi(V) = Ax(V) forapair (i,k) with i Æ k 


would have to exist because of continuity reasons. But this would contradict the strict 
hyperbolicity on N. 


The eigenvectors s;(V) (i = 1,2,- -- , m) belonging to the eigenvalues \;(V) (i = 
1,2,- -- , m), respectively, are linearly independent. 


Let V (x,t) be a solution of (2.1) which is smooth on a domain G C Q and with 
Vz, t) EN , V(a,t) EG. 


The i-characteristics x = zalt) corresponding to this solution (cf. (1.18)) are 
uniquely determined due to the Picard-Lindelgf theorem if the functions \;(V(z, t)) 
fulfill a Lipschitz condition with respect to x (i = 1,--- ,m) and as far as initial 
States £(;) (to) with (a (;)(to) , to) € G are prescribed. Of course, there are also exis- 
tence and uniqueness theorems with weaker assumptions available, e.g. the existence 
theorem of Peano combined with the uniqueness condition of Nagumo. 
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Then, along an 7-characteristic, 


d 
gV to: t) 

= — [Jf (V (a(t), t)) — A (V (z(t), t)) - T] V (ealt) t) 63D 
where T in this context means the m x m unit matrix. 


In case of m = 1, the righthand side of (3.31) vanishes as already observed earlier, so 
that V = v is constant along a characteristic which therefore becomes a straight line. 


Does an analogon exist for m > 1 ? In other words, does a weak solution V and 
does a set of i-characteristics z; belonging to this solution exist so that the solution 
is constant along these characteristics which, thus, become straight lines? 


A positive answer requires in case of ôs V (x; (t), t) # 0 the validity of the relation 
3V (zalt) t) = si(V (zt) t) VtE>0”. 


Solutions of this type can really be specified: If we introduce so-called centered 
waves V = V*, centered at a position Po = (x0, to) € G, i.e. waves of the form 





V (x,t) = V* G — z) (3.32) 
t — to 


with 

V*(& E€ NAO[CI(R]” VP=(2,t)eG , é:=—— 
These centered waves solve (2.1) if they fulfill 

OV + Oe f(V) = (rE J f(V*) + HET) V” (E) = 0 
which can also be described as 


(J FOE) —€ 1) VE) =O. (3.33) 
Hence, for one i € (1, 2,--- , m), the pair of relations 


E= ài (V*(6)) (for all possible values of £) (3.34) 


must hold. 


>*Here, the norm of the vectorfield V is chosen suitably. 


86 3 Entropy Conditions 


In this case, differentiation of the second equation of (3.34) with respect to € together 
with the first equation yields (cf. footnote 51) 


= (WvAi(V), si(V)) . (3.35) 


Definition: 
If (VyAi(V), si(V)) # 0 holds for all V € N witha set N C R”, the i-th 
characteristic field (A;(V) , s;(V)) is called genuinely nonlinear on N. 


Remark: 
Because the scalar case implies A(V) = X(v) = f'(v) which yields 


VyA(V) = f"), 


(3.35) can only be fulfilled if f"(v) # 0, and s(v) in this case is a scalar factor 
s(v) different from zero. It can be normalized by s(v) = Fo which makes (3.35) 
equivalent to f’(v) > 0, i.e. to strict convexity of the flux (or, if we normalize by 
TONE to strict concavity). 


We are therefore going to look upon (3.35) as the transition announced earlier of the 
strict convexity of the flux as defined for scalar problems to the i-th characteristic field 
incase of m > 1. 


Let us now assume N to be a simply connected domain in R™. Moreover, let 
V°? €e N (V° 40) bea given state and ĉo := \;(V°) . If & will be used as 
starting point with initial condition 


V*(o) = V°’, (3.36) 
the second equation in (3.34) will obviously be fulfilled at this point £o. 


The first equation of (3.34) represents a system of ordinary differential equations 
where (3.36) is used as its initial condition. We know that there is a unique solution of 
this system in a certain neighbourhood of (ĉo, V°) provided that the assumptions of 
a standard theorem like the Picard-Lindelgf theorem are fulfilled. This solution runs 
with respect to V in N, let’s say for éo —-a < € < éo +a with sufficiently small 
a > 0. Also (3.35) then holds for V = V*(&) because of the assumptions. 

Since also the first equation of (3.34) was fulfilled, also (3.33) holds. Integration of 
(3.33) together with (3.36) then yields the fact that the second equation of (3.34) is 
fulfilled, too, at least for the values of € under consideration and if a will be again 
somewhat diminished if necessary. 
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We conclude that there is a nonvanishing smooth solution V, namely a wave centered 
at Po, belonging to a given nonvanishing smooth initial state V° and to a particular 
integer 7, leading to an 7-th characteristic field which is genuinely nonlinear on N: 














V (x,t) = V* k — z) o V*(A(V®)) = V° (3.37) 
— to 
for all (x,t) € Q with 
\(V°) -a< — <i(V°) +a with a>0 sufficiently small. 
— to 
(3.38) 
Along each of the straight lines 
T — To 
— = c= const 
Paty 
with 
ri(V°) —-a<c< r(V°) +a 
forming a fan, the solution 
Vv" E — z) = V*(c) 
t — to 
is constant. 
t 
om x0) 
N 
W 
AD 
oe (x—X0) 
Po= (xo, to) _ 
x 


Figure 17: Fan of i-characteristics of centered waves 
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Thus, along these straight lines given by 
za (t) = zo + c(t — to) , 
the relation 


d do 
LV o0) = FV") =0 3.39) 


holds. But because V*’(c) is an eigenvector belonging to \;(V*(c)), this also applies 
to 


1 x 
On V (x(t), t) = rere He) x 


It follows from (1.18) and (3.33) that the straight lines generating the fan and belong- 
ing to this particular integer 7 are the 2-characteristics formed by means of the solution 
(3.32) through the point Po, and the centered waves constructed so far are constant 
along these lines. 


Let us now assume that WV! is a state in the neighbourhood of V°. Then, by 
continuity arguments, \;(V") lies in a neighbourhood of \;(V°) such that 


IAV?) = AV’) <a 
can be fulfilled for all states V! sufficiently close to V®. 
If particularly 
NOV) < Ai V!) < AVO +a? 
we find the result 
v* k = z) _y? 
toto 


along the straight lines 





pai cesu) 
— to = —— r (T — 20), 

me NOVO A 
so that the solution V*(€) is constant along the lines and equals V° if this is its 
value at ĉo. The same arguments show that the solution equals V! along the straight 
lines 


if this is its value prescribed at €. 


In other words: In case of A;(V°) < 2,;(V'), the state V° can be smoothly 


connected from the right with the state V! by a centered i-rarefaction wave™. 


In case of As(V") < Ai(V°), change the roles of V+ und V°. 
“of. also the definition that follows formula (3.43) 
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Remark: 
In the scalar case and if f is strictly convex, \;(V°) < \;(V') means f'(v?) < 
f' (vt), hence 


v? < vt. 


In other words: If t is fixed, the states decrease for decreasing x, i.e. from the 
right to the left (rarefaction). Indeed, if the entropy condition (3.20) is fulfilled, this 
rarefaction can only be caused by a continuous process because of vg = vy. 


Assume now the system (2.1) to have discontinuous solutions besides or instead of 
smooth rarefaction waves, and let us assume that these discontinuities form certain 
curves I’ along which the Rankine-Hugoniot condition is satisfied. 


We ask again for conditions easy to apply in order to characterize uniquely the phys- 
ically relevant entropy solution”®. 


Here, we restrict ourselves to problems whose characteristic fields are genuinely non- 
linear on N = R”, 


A weak solution will be called an entropy solution if it respects a condition of the 
following type: 


— If initial values are prescribed along a non-characteristic curve, the condition will 
only be fulfilled by one weak solution. 


— The condition coincides with (3.24) as far as scalar problems are concerned. 


If [V] = V~—V,. are the jumps along a curve [ of discontinuities of a weak 
solution moving along the x-axis with speed @, let k be the particular integer with 
Ak(Vr) < Ô< Ak41(Vr) , KE {1,2,---,m—-1}, (3.40) 


or ô< (V) or ô > Am(V_)*. 


In order to determine the m components of V, at Po € I from the given initial 
values, only the values of Vg along the particular characteristics through Po are 
available which show the steeper ascents 


1 1 


1 
ee ees ts Se (cf Fig. 18)°? j 
An(Wo) | An—1(Vo) A1(Vo) Pere 


STf the flux f is a concave function like the flux of the traffic flow model, the word ’left’ has to be 
replaced by ’right’ and ’decrease’ by ’increase’. 

>°Remember that condition (3.12) could not necessarily be applied in case of more than two equations. 

>7analogously to scalar problems, i.e to m = 1 where only convex fluxes were taken into account 

Sof. (3.30) 
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If exemplarily a system of decoupled equations is under consideration®, along each 
of these k steeper characteristics one and only one of the m components of V (Po) 
can be determined. Hence, m — k informations are missing. 


If simultaneously 
Aj (Ve) <d0< Aj+1(Ve) y 


holds, corresponding arguments lead to the result that there are only m — j in- 
formations available for the computation of the m components of Vy. Hence, j 
informations are missing. 


t 





corresponds to À% AKA x 


Figure 18: Sketch concerning the completeness of informations needed in order to 
ensure uniqueness 


Thus, m — k + j informations are missing for the unique determination of the two 
quantities V, and V,. But m additional informations concerning the relations 
between Vy», Vy and 0 follow from the Rankine-Hugoniot condition (2.15). One 
of them must be used for the elimination of ô. Thus, the real number of missing 
informations is only 1 — k + j. In other words, uniqueness can only be expected for 


k=jr+l, 
i.e., if the relations 
MVr) <6 <Angi(Wr) >  An-1(We) <ô < An(Ve) (3.41) 


hold simultaneously. 


»The particular choice 0 € T in Fig. 18 does not lead to a loss of generality. 
cf. the linear situation (1.19) 
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Although this result was only illustrated by an exemplary situation, we accept it as a 
hint on a general definition of an entropy solution: 


Definition: 

V is called an entropy solution if there is an index k for which (3.41) holds. A 
discontinuity of this type is then called a k-shock and the inequalities (3.41) are called 
entropy inequalities or (Lax-) shock conditions. 


(3.41) includes 
Ak(V 7) << Ak(V o) ; 


so that (3.41) really coincides with (3.24) as far as scalar problems are concerned. 
Thus, both properties claimed earlier with respect to the formulation of an entropy 
condition for weak solutions of systems of conservation laws are respected. 


There are some connections between an i-characteristic field and the so-called i- 
Riemann-invariant: 


Definition: 
A scalar function w = w(V) defined for all V € N C R” is called an i-Riemann 
invariant on N, if 


(Vyw(V), sV) =0 YVEN. (3.42) 
Remark: 
Thus, the property of an %-characteristic field to be genuinely nonlinear, coincides 


with the property of 4;(V) not to be an i-Riemann invariant. 


If the graph of a centred i-rarefaction wave V is situated in WN, its 7-Riemann- 
invariants are constant. This follows from (3.34): 


0, w (V(x, t)) = (Vyw(V) , 8i(V)) yav«(a,) =0 (3.43) 
and similarly Oj w (V (x,t)) = 0. 
Definition: 
If a solution V of the system (2.1) belonging to any initial function Vg is smooth 


on a region G C Q., and if all the i-Riemann-invariants w(V (x,t)) are constant on 
G, V is called an i-simple wave or an i-rarefaction wave on G. 
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3.4 The Ansatz due to Kruzkov 


The way used by S.N. Kruzkov®! in order to introduce the idea of a weak solution 
of problem (2.1) seems to differ from the idea of Lax. In order to understand the 
situation, it is sufficient to restrict the presentation of Kruzkov’s Ansatz to only one- 
dimensional scalar problems. 


Also Kruzkov tries to find his weak solutions in the space L° (Q) as it was already 
done in the Lax definition of a weak solution. Moreover, also the test functions used in 
(2.2.) by Lax play a crucial role in Kruzkov’s definition, too, with the one restriction 
that these functions now vanish along the whole boundary of their compact support, 
i.e. also along parts of the boundary belonging to the x-axis: ®(x, 0) = 0. This leads 
to the necessity to pay regard to the initial function in another way than by the second 
integral in formula (2.2). 


Kruzkov calls a function v € L° (Q) a weak solution of the problem 
Op + Oxf (v) =0, v(z, 0) = voz) ) 


if the following two conditions are fulfilled: 


/ faatt) E E 
Q 


OE E sol} d(2,t) >0 844 


y € CEQ) with © > 0 and with (7,0) =0, Yc € R, 
and 


R 
lim |u(a,t) — uo(x)| dx = 0 VRER (teb, TS (3.45) 


t—0 -R 


for a value of T > 0 and for a set € of measure zero so that v(x, t) is almost 
everywhere well defined as a function of x for every fixed t € [0,7]\¢é. 


Kruzkov could show that there is a unique weak solution in L'*(Q) in the sense of 
(3.44), (3.45) provided that the flux f is smooth and strictly convex. Smoothness and 
strict convexity of the flux were assumptions also made by Oleinik and Lax. 





6! Soviet Math. Dokl. 10 (1969), 785-788 
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Hence, in Kruzkov’s definition, the characterization of an entropy solution does not 
consist of an equation (2.2) to be solved for all test functions together with an addi- 
tional entropy condition (3.12) but only of an inequality®. 


One can say that Kruzkov’s definition already includes an entropy condition. 


Indeed, as far as the scalar problem is concerned, the left side of (3.44) corresponds 
with (3.12) if ®(2,0) = 0 is taken into account, if for every particular constant 
c E R the entropy functional (v) will be identified with |v — c| and the entropy 
flux F(v) with sgn (v —c) [f(v) — f(o). 


If S and F are chosen in this way, also the connection between these two functions 
demanded by (3.13) is kept piecewise, namely for all values of v with the exception 
of v = c. Moreover, the particular functional S is also convex as entropy functionals 
were expected to be but no longer strictly convex, unfortunately. 


On the other hand, there is a somewhat stronger demand than formulated by Lax, 
namely that the validity of (3.44) has not only to be ensured for all nonnegative test 
functions ® but at the same time also for all c € R. Thus, (3.45) must now hold for 
all test elements 


de {(O,c) | GEC), B>0, G(w,0) =0; ce R}. (3.46) 


It can be shown by more than one chain of evidence that Kruzkov’s weak solution 
coincides with the Lax-Oleinik entropy solution provided that the flux is smooth and 
strictly convex. One of these chains follows from numerical aspects: There are finite 
difference approximations whose numerical solutions converge for decreasing step 
sizes to the Lax-Oleinik entropy solution but also —as can be proved by a slightly 
different convergence argument— to the Kruzkov solution®. 

It should be noted that there are problems from fields of applications where the as- 
sumptions on smoothness and strict convexity of the flux are violated. But Kruzkov 
together with Panov could show that a unique solution of the problem (3.44), (3.45) 
does even existif f is merely continuous™. Thus, in the scalar situation, the Kruzkov 
concept of an entropy solution is more far-reaching. 


But one has to pay for the renunciation of convexity and smoothness: It can happen 
that one of the properties of the entropy solution vg important from the point of 
view of numerical methods gets lost, namely that an initial function vg with compact 


®and of a condition that respects the initial condition in a suitable way 

& A non-numerical argument concerning the connection between both definitions can be found in the 
paper already mentioned in footnote 38. 

“S.N. Kruzkov und E.Y. Panov: Soviet Math.Dokl. 42 (1991), 316-321 
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support leads to a compact support of vg(-,t) (for every fixed t > 0) as it is the case 
if f is smooth and convex. A test example for this loss is realized by 


Q 
fœ) = KE with 0<a<l, 
a 
(3.47) 
0 for x <-—l1 
volz) = 1 for -l<2x<0 
0 for xr>0 
The solution reads 
0 for t>a(x+1) 
1 fr xz<t<alxr+1 
vp(x,t) = a ry (3.48) 
(<) for t<a 
x 
(cf. Fig. 19). 








Figure 19: Kruzkov-Panov solution (3.48) for 0 < t < -= 


l-a 


Whereas the smooth and strictly convex case behaves hyperbolically, a merely contin- 
uous flux can lead to a more parabolic character of the problem. Numerical procedures 
for problems of this particular type were under some suitable conditions recently be 
developed for the first time by M. Breuss®. Here already the Kruzkov-Panov exam- 
ple (3.47) shows that the important CFL-condition (cf. sections 6.1, 7.3) can not be 


6M. Breuss: Numerik von Erhaltungsgleichungen in Nicht-Standard-Situationen, PhD-Thesis, Dept. 
of Mathematics, University of Hamburg, 2001 
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fulfilled as far as explicit finite difference methods (cf. sections 4.1, 6.4, 7.2) are used. 
Thus, implicit methods have to be taken into account, and this can really be done 
successfully. 


Examples of mathematical models of real world problems in fluid mechanics where 
standard assumptions like strict hyperbolicity, strict convexity or smootheness of the 
fluxes are not necessarily fulfilled occur in the theory of flows through porous media, 
e.g. if oil production is concerned: 


The scalar so-called non-standard Buckley-Leverett equation 


2 


Bw + ð ( ) =o,aer 


v 
v2 +a(l — v)? 


models the displacement of an oleic phase by an aqueous phase in a porous medium. 
In this situation, the flux is neither convex nor concave. 


Tveito and Winther® considered theoretically with respect to stability the system 


Os + Onf(s,c) = 0 
O.(sc+a(c)) + Or(cf(s,c)) = 0, 


which also models the displacement of oil by water in a porous medium. But in this 
case, the water is thickened by dissolved polymer in order to increase its viscosity. 
The difficulty of this example arises from the fact that the eigenvalues of the Jacobian 
are real but can coincide so that the system is hyperbolic but not necessarily strictly 
hyperbolic. 


Proc. 3rd Intern. Conf. Hyperb. Probl. (B. Engquist, B. Gustafsson, ed.), (Chartwell-Bratt, Upp- 
sala), 888-898 (1991) 


4 The Riemann Problem 


4.1 Numerical Importance of the Riemann Problem 


Treating the traffic jam dissolution example we met the first time with a Riemann 
problem. It consisted of a conservation law of type (2.1) together with an initial 
function being piecewise constant with certain jumps at isolated positions. 


One tool constructed in order to solve differential equations numerically is the so- 
called Finite Difference Method (FDM) where one tries the unknown values of the 
exact solution at isolated points to be approximated (cf. also chapters 6, 7). 


Particularly in case of problem (2.1), we ask for numerical values vi at isolated 


positions (x), tn) of the upper half plane Q expected to be good approximations for 
the values V g(2,,t,,) of the exact entropy solution at these points. 


A particular possibility for the construction of such FDMs consists of covering the 
positive time axis in a first step by a not necessarily equidistant grid {t, | (n = 
0,1,---)} of time step sizes Atn := tn+1 — tn together with fixing isolated grid 
points x, ; (v =0,+1,+2,---) along the x-axis. Thus, by the cartesian product of 
these two sets, 2 will be covered with a net of grid points parallel to the axes. Also 
the step sizes hy := £y+1 — x, don’t need necessarily to be equidistant (cf. Fig. 20). 














If the computation of the values yor along the time level t = t,41 does only 


use the knowledge of the values vir along the time level before and computed in 
the very last step, we call the scheme a one step method. Multi step methods also 
include approximate values along earlier levels, but we are going to restrict the present 


survey to one step methods. 


Let to = 0 be the level to start from with approximate initial values 


5 (ay +av41) 
/ Vo(a)dx (v=0,+1,+2,---). (4.1) 


4 (ty —1+@v) 





2 
V9) = ——=___ 
Ly4+1 — Lyp-1 
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Figure 20: Grid of an FDM 


It seems to be reasonable to describe the difference between the approximate solution 
and the unknown exact solution in terms of the topology of the space forming the 
basis of the original problem, which is in the particular situation under consideration 
the space L loc (gy, This difference or truncation error should be estimated using a 
suitable norm. Hopefully, these errors will decrease for decreasing step sizes and will 
tend to zero sufficiently fast if the step sizes tend to zero. If this holds, the method is 
called convergent. 


Obviously, a direct comparison of the approximate solution to the exact solution in 
terms of the toplology of the original space can take place if the discrete values 


(n) 


ò \ of the approximate solution defined only at isolated points can be extended 
to the inter grid points suitably, e.g. by 


V (x,t) = v) for a <a< Titoki, tn < t< tni; 


V 


(4.2) 





(v = 0,1, +2,--- ; n =0,1,2,---). 











By this special extension or reconstruction the discrete function now becomes a func- 
tion V constant on each of the described rectangles. Hence, V isa step function 
on Q and, thus, an element of loca), 


If we restrict V to the time level t = tn, this restricted function V (a, tn) of the 
space variable x is a step function, too, from which the one step method together with 
the reconstruction now generates V(x, tn41). 


If V (a, tn) is used as an initial function for the original problem in order to generate 
a weak solution of the differential equation in (2.1.) for the region t > tn, a Riemann 
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problem comes about. Its exact solution if particularly restricted to the time level 
t = tn41 would lead at the grid points along this level to values which could be 
looked upon as approximations to the weak solution of the original problem at these 
new level grid points. This idea was for the first time realized by S.K. Godunov’ and 
was later often used by other authors as a basic idea for further numerical schemes. 
Herewith, effective approximate Riemann solvers were developed including local 
linearizations. 


We are therefore now going to study the Riemann problem at least in case of linear 
problems with constant coefficients. 


4.2 The Riemann Problem in the Case of Linear Systems 


Let us recall problem (1.19), i.e. 
OV + AOdO,V =0 F 
with a constant (m,m)-matrix A, and let us prescribe initial values 


= V; for «<0 
Vole) = { V, for «>0 oe 


with constant vectors Vy, und V, different from each other. It causes no loss of 
generality using the origin as the position of the jump. 


Assume the problem to be again strictly hyperbolic so that the eigenvalues A; of the 
matrix A are real and different from each other. The eigenvectors belonging to these 
eigenvalues are then obviously real, too, and linearly independent. Thus, they can be 
used as a basis in JR”, and in areas were the solution is smooth, (1.21) and (1.22) 
yield 


V(2,t) = X big(a — À; t) si (4.4) 
i=1 


with 
Volz) = (io (2), b29(2),+++ , mo (#))" = S™'Vo(2) . 
This leads particularly to 


Vo (x) =; (w1, w2, aen Wm)? = SiV, for «<0 
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and analogously to 
Vo, (x) =: (w1, wa, was Wm)? = StV, for 2 >0 
with certain constants w; and w;. 


If the eigenvalues are ordered as in (3.30), 


k m 
V(z,t) = `i Si t 5 Wi Si (4.5) 
i=1 


i=ko+1 


follows from (4.4) provided that for a given position (x,t) the first kg eigenvalues are 
the particular ones for which x — A;t > 0 (i = 1,2,--- , ko), whereas x — Ajt < 
0 (¢=ko+1,ko +2,--- ,m) holds for the other m — ko eigenvalues. 


Let us now look at the 7-characteristics x = A;t arising from the jump position and 
let to > 0 be fixed. Let us move for this fixed tọ from the left to the right along 
a line parallel to the x-axis. This line crosses the ko-characteristic at the position 
Xo = Ako to (cf. Fig. 21). As soon as this will happen the coefficient ko in formula 
(4.4) has to be replaced by wko- 





Figure 21: How to solve the linear Riemann problem 


The solution V described by (4.5) then obviously jumps. The height of this jump is 


[V] ko (T0, to) = (Wko — Ūko) Sko (4.6) 


and leads to 


A [V] ko (xo, to) = (Wko z Who) Ako Sko = Ako [V]ko (£0, to). (4.7) 
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x—-A3t x—-Ant x—Àıt 


Figure 22: Characteristic cones of the solution of the linear Riemann problem 


We know already that discontinuities move along characteristics as far as linear prob- 
lems are concerned. This makes the ko-characteristic line through the jump position 
a ko-shock Iko if Wko Æ Wk so that (4.7) just becomes the Rankine-Hugoniot 
condition (2.15) applied to the special linear situation. 


(4.5) and (4.7) yield the result that the solution of our Riemann problem for t > 0 
can also be formulated as 


V(2,t)=Ve- X` (wi — ùi) si (4.8) 
y2 
or as 
V (x,t) = V,+ bp (wi — wy) Si (4.9) 
ysg 


(cf. Fig. 22). 
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5.1 The Navier-Stokes Equations Model 


An ideal fluid was defined to be a material without friction between neighbouring 
particles of different flow velocities, hence without tangential forces along the surfaces 
of arbitrary volumes W(t) picked out of the fluid flow. The mathematical model of 
ideal fluids can therefore only describe certain particular situations, e.g. flows of 
low speeds, low densities, low viscosities or high temperatures etc., and its results can 
seriously contradict the reality if such properties are not ensured. The forces acting on 
bridge piers parallel to the flow of calmly flowing small rivers were already mentioned 
when the Kutta-Zhukovsky formulas were presented in section 1.3. 


Thus, if the tangential forces generated by the viscosity play an important role in a 
real flow situation they can’t be neglected in a mathematical model expected to be 
able to describe the flow in a satisfactory manner. 


Let us therefore try to revise the Euler equations model of ideal fluids by taking the 
viscosity forces additionally into account in order to generate a model of real fluids. 


Reminding the reader to the fact that viscosity in our context always means friction 
between fluid particles®’, these additional forces can only occur in case of motion, i.e. 
for non-vanishing velocity of the flow, and lead then to an additional term within the 
conservation law of momentum: 


The righthand side of (1.5) consisted of the forces per unit volume k only, and the 
part of these forces falling to the volume W (t) therefore was given by 


l kica: 


W(t) 





of. (1.54) and its explanation 
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Forces generated by friction of the fluid particles on the outside of this volume but 
moving along its surface W (t) are obviously proportional to the area of W (t) 
and have tangential directions, whereas the forces generated by pressure of the outside 
fluid are proportional to this area, too, but directed perpendicular to OW (t). 


If s = (s1, 82,83)! are the forces per unit area generated by the outer fluid and acting 
on the surface OW (t) of an inner volume W(t) arbitrarily picked out of the fluid, 
the relation 


s = s(t,x,y,z n), (2,y,z) € OW(t) (5.1) 


holds where n is the unit vector normal on OW (t) at (x,y,z) € OW(t) with 
outward direction. 


Herewith, s depends linearly on n: 
s = (oij) n (5.2) 
with oij : (x,y, Z, t) > o4j (x,y, z, t). 


(5.2) follows from the fact that within the conservation law of momentum, now to be 
written as 


i faat (Zav) q + div (<2) a} d(x, y, z) 


W(t) 
= if kd(az,y,z) + | sdo , 


W(t) aw (t) 
the volume terms vanish of a higher order than the surface terms in case of W(t) — 0. 
Here, the stress tensor (c;;) is treated like a real (3, 3)—matrix. 
Moreover, Euler’s law of vanishing momentum, i.e. 
J tr mldey2) + ff ry sdo = 0, (5.3) 
W(t) OW (t) 
r = (x,y, z)!, yields the symmetry of the stress tensor: 
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Vanishing friction in case of ideal fluids has to be modelled by o;; = 0 for i Æ j. 
Taking (5.2) into account, this leads to 





O11 022 033 : —p (pressure) 





because the vectors s and n are parallel in this situation. This is a slightly more 
general explanation of (1.5.). 


As far as real fluids are under consideration, it can be expected the terms oj; (i 4 j) 
to be small, and the terms oj; will certainly differ from —p only to a small degree. 
We therefore put 


a. Tii tp for i= 7 

Tij = i A : (5 5) 
so that all the o;; are small. 
Because of the symmetry of (0;;), also the tensor (G;;) is symmetric. 
If the exterior forces k per unit volume are replaced by the forces k per unit mass, 


/ {aa + (žav) q + div (<2) a— pk} d(x,y, z) 


W(t) 


p e I pnis f (Gy,) do 


aw (t) aw (t) 
follows. 


Obviously, if a fluid does not move or if neighbouring particles do not move with 
different velocities, friction can’t occur. In modelling the flow, we therefore assume 
the terms 6;; (i,j = 1,2,3) to depend only on the components of the vectors 
V u; (i = 1, 2,3). Here, we restrict ourselves to so-called Newtonian Fluids, i.e. the 
terms ij depend on the components V u; homogeneously and linearly: 


3 
bp = SO) Only (5.5) 
pv=l 


Thus, every oj; is a function of the elements of the Jacobian Jw: 
Õij : Ju — õij(Ju) 5 (5.6) 


and the ali) are constants. 
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Analogously to the case of ideal fluids where every axis is a main stress axis, we 
also expect in our model the relation (5.6) to be independent of orthogonal linear 
transformations, i.e. 


Al (6i;) A= (aij(ATJuA)), VAER®) with ATA=T. (5.7) 
It is trivial that the matrix B := divu-J is symmetric, and because it is even a 
diagonal matrix whose diagonal elements equal eachother it shows the property 
A'BA=B. 


The elements b;; of B depend on the elements of the Jacobian, i.e. B = (bj;(Ju)). 
Because of div u = trace Ju and because traces of matrices are invariants of orthog- 
onal transformations, also 


(bij (A7 Ju A)) = (by(Ju)) = B 

follows. 

Thus, B fulfills all the expected demands on (6;;). But also 
D := Ju + (Ju) =: D (Ju) 


is symmetric, depends on the elements of the Jacobian in a homogeneous and linear 
way and shows the property 


ATDA = A" JuA + A" (Ju) A 





= A" JuA + (AT JuA)', 
i.e. A'D(Ju)A = D((A'JuA)) . 


Thus, also D fulfills all the demands on (6;;), and therefore also AB + nD with 
arbitrary constants A,n does so. It can be shown that these are already all types of 
matrices which satisfy our demands. Hence, 


I {aa+ (247) q + div (54) a- pk} d(x, y, z) 
W(t) 
(5.8) 


— / (-pl+rAB+n7D)ndo=0. 
aw (t) 
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By means of the divergence theorem, 


pndo = Í V p d(x, y, z) 
aw (t) W(t) 
div b; 
and f Bndo = div bə | d(z,y,z) 
aw (t) W(t) div bs 


follow where the b; are the row vectors of B, i.e. 
Bndo = I V (div u) d(x, y, Z) , 
aW (t) W(t) 


and similarly 


div dı 
I Dndo = / divd, | d(x,y,z) 
aW (t) Wit) div d3 
where 
divd; = 20, (Oru1) + Oy (yur + ruz) + Oz (Ozu1 + ðru3z) 
= Auw + rru + Oxy tla T OnzU3 








= Au + (divu). 


Analogous formulas hold for div dz and div d3. 


The equation (5.8) therefore reads as 
1 : 1 s 
/ [aas (20.7) q + div (<a) q — pk 
P P 
W(t) 


+ Vp — rAVidivu) — nAw — nian) bale = 0 


for every part W(t) picked arbitrarily out of the volume of the fluid, i.e. 


O,(pu)+(u, V) (pu) +(pdiv u) u—pk+V p—(A+n) V(divu)—n Au = 0 
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or 
pu + put lu, V p) ut+p(u,V) ut (pdivu) u—pk+Vp 


—(A+n)V(divu) -nAu = 0. 


Because of pdivu + (u, V p) = div(p u), the relation 





{Op + div(pu)} utp utp lu, V) u—pk+V p—(A+n) V (div u)—7 A u = 0 


follows, and taking the continuity equation into account, we conclude after division 
by p with the so-called Navier® - Stokes”? equations 





À 1 5 
ðu + (u, V) u — T Vidiv) -— nau =- Vp +k]. 69) 





Here 7 is called the first viscosity coefficient. Of course, historically there were 
several types of arguments leading to the Navier - Stokes equations, and this led to the 
convention not to call A the second viscosity coefficient but € := À + Zn. These 
coefficients are constants depending on the fluid material as well as on its temperature. 
The unit of measurement is g- cm! -sec~! = Poise (cf. footnote 71). 


In case of incompressible flows characterized by divu = 0, A does not occur in 
the equations, and if also the exterior forces k do not play a role, the Navier-Stokes 
equations are reduced to 


1 
gu + (u,V)u-vAu+--Vp = 0, (5.10) 
p 


where v = 3 is called the kinematic viscosity. 


By similar arguments, the law describing the conservation of energy takes for viscous 
flows the form 


OE + div ((E + p) u) — A - div(div u - u) — n- div (Ju + (Ju)") u) = 0. 
Some of the terms within this equation can be simplified, e.g. 
div(divu-u) = (V(divu),u) + (div u)? 


etc. 


Claude Louis Marie Henri Navier (1785 - 1836); Dijon, Paris 
George Gabriel Stokes (1819 - 1903); Cambridge 
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The Navier-Stokes equations model is excellently confirmed by experiments. 


An impression of the magnitude of the kinematic viscosity and its sensibility with 
respect to temperature is given by the following table: 


0.00125 
0.00123 
0.00117 


Machine oil 
(depending on the 
brand; approx.) 





In order to study the behaviour of a real ship in the ocean or of a wing in air flow etc. 
by means of experiments, engineers normally use a small model of the ship dipped 
into a wave canal or a small model of the wing brought into a wind tunnel, respectively. 
Let us ask for conditions under which the results of the experiments are in accordance 
with the reality. To answer this question it is advisable to use dimensionless quantities. 


Let us treat for convenience the case of an incompressible flow: 


Assume up = (ug,0, 0)? to be the flow velocity in a great distance of the solid, i.e. 
of the hull of the ship, of the wing etc., let L be a characteristic length of the solid, 
e.g. the length of the ship, the length of the wing’s airfoil etc., let T' be the time the 
fluid needs to cover the stretch of length L with velocity uo, and let P = 5 uo? be 
the dynamic pressure in a great distance of the solid. 


Let us consider the particular case where exterior forces per unit volume can be ne- 
glected and let us now introduce the dimensionless quantities 
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The incompressible Navier-Stokes equations then turn out to become 


y 7 i 
oy + (1,7) a+ Vp - = Aa = 0, (5.11) 
where 
L L 
Rea ee (5.12) 
n V 


is the so-called Reynolds number and 


V := (0%,05,0s)" , A = Oza + Ogg + Oxz 


Thus, the experiment coincides with real world behaviour of the solid put into the flow 
if the Reynolds numbers of both situations coincide. This leaves a great freedom with 
respect to the particular choices of uo, L, v realized in the experiment. The flows 
are then called similar. 


In case of n Æ 0, the Laplace operator occurs in (5.10) so that this system of equations 
now becomes a parabolic system if instationary flows are concerned or an elliptic 
system if higher dimensional stationary flow is under consideration, respectively. 


The continuity equation as well as the energy conservation law already known from 
the model of ideal fluid flows remain unaffected by the viscosity and remain valid. 
This also holds for the additional equation of state as far as gas flows are concerned. 


There are only a few situations where the Euler equations or the Navier-Stokes equa- 
tions can explicitly be solved by wellknown functions. In all the other situations one 
has to solve them approximately by means of numerical procedures. 


Two of the problems which can be treated explicitly will be studied in the next section, 
and both problems are examples of stationary flows. 


Remark: 

The Navier-Stokes equations do not take the particular additional buoyancy phenom- 
enons into account arising in compressible gas flows from density variations caused 
by temperature differences between different parts of the gas. In order to respect these 
phenomenons, the Navier-Stokes equations have to be enriched by additional terms. If 
the enriched equations are also written in a dimensionless form, further dimensionless 
parameters besides the Reynolds number occur, particularly the Prandtl number 


v 
Pr := \ PI 


with the specifiic heat conductivity A. 
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5.2 Drag Force and the Hagen-Poiseuille Law 


Viscous flows can be laminar or turbulent. It is laminar if neighbouring fluid par- 
ticles move more or less parallel to eachother along stable trajectories, whereas tur- 
bulent flows are characterized be a disordered flow that superimposes the main flow. 
Thus, the turbulent flow is an instability phenomenon. When Reynolds studied the 
flow of liquids through small glass tubes, he noticed that the transition from the lam- 
inar to the turbulent state of the flow takes place in a very abrupt way depending on 
the parameter that is now called the Reynolds number (cf. (5.11)). 


The transition to instabilities will briefly be discussed in the fifth section of this chap- 
ter. In this section the flows are assumed to be laminar. 


Our first simple application of the Navier-Stokes equations concerns the laminar vis- 
cous stationary flow in y-direction of small velocity parallel to a given area ôF of 
a plain solid plate where this area is part of the (y, z)-plane. Let the fluid be incom- 
pressible in the sense of p = const. Because of the friction between the fluid particles 
moving along and the fluid particles sticking to the rigid body in accordance with the 
no-slip condition, there is a force acting on ôF. Obviously, this force has only one 
component in y-direction, and the velocity is independent of the space variables y and 
2: 


u = (0, u(x), 0)" with u(0)=0 (no-slip condition). 


The undisturbed velocity of the fluid that passes the surface area at a great distance is 
assumed to be constant, i.e. O,u(z) - 0 for z > o. 





Figure 23: Flow along a plain membrane 
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Compared with (5.10), the forces per mass unit occuring within the Euler equations 
are now enriched by the additional term v A u so that the forces per volume unit are 
now modelled by 


0 
w= } Open) |. (5.13) 
0 


The resistance of the relevant area acting against the flow in a direction parallel to 
itself and arising from the friction or —vice versa— the drag force acting on the plate 
by the friction of the fluid comes therefore up to be 


ee 0 
oP | war = —nuz(0) F 
0 0 


Hence, the drag related to the friction of a laminar flow is given by 


W = n 0,u(0) OF. (5.14) 


Coming back to Reynolds’ experiments, let us study as an example the one-dimensi- 
onal laminar stationary cylinder-symmetric flow of an incompressible liquid through 
a circle-cylindrical pipe of radius R where the gravity force is neglected, where other 
outer forces do not occur and where the no-slip condition u = 0 holds at the interior 
surface of the pipe (cf. Fig. 24). 





Figure 24: Laminar pipe flow 


Because of u = (u, 0,0)" the relations 


(u, V) u = (udu, 0, 0)” 
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and 
Au=(Au, 0,0)", divu = ôu 
follow. 


Taking the continuity equation div u = 0 into account, we find 0,u = 0, thus also 
zzu = 0. Hence, (5.8) reduces to 


Au 1 ozp 0 
stiro khta ap E a (5.15) 
P\ 0 P \ ap 0 


This leads to 


yp = ðzp=0, ie. p= p(z), 
and therefore to 
O 
Au = ðyyu + ðu = ZE, (5.16) 
n 


The left side of (5.16) is represented by a function that depends only on the variables 
y and z, whereas the right side depends only on x. So, 0,p has to be constant. 


The particular case of a flow considered here can therefore only occur if p depends 
linearly on z. 


If polar coordinates (r, p) are introduced in the (y, z)-plane, and if the cylinder sym- 
metry is respected, i.e. u = u(r) independent of y, the well-known transformation 
rule 


R E E Oa 
r 


re 
transforms (5.16) into 


ð, (r ðu) = Ôe ,, 


Because p is independent of (y, z), hence independent of (r, p), integration leads to 


O,p r? 
Oi a Ae (5.17) 
n 2 


with an integration constant cy. 
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Besides the no-slip condition u() = 0, also the boundary condition 
(A-4)-<0] < œ% 
has to be fulfilled so that cı has to vanish. 


If (5.17) is now integrated, 


Oxp r’ 
ur) == = +c 
NF ote 
results where the constant c2 has to be chosen in such a way that the no-slip condition 


is fulfilled. This finally yields the so-called Hagen-Poiseuille law 
Orp 2 2 
= —-— (R — ; 5.18 
OES ga ete) (5.18) 
For every fixed x, this law shows a parabolic velocity distribution, and the flow follows 
the direction of decreasing pressure (cf. Fig 25). 


v (r) 4 





























Figure 25: The Hagen-Poiseuille flow through a pipe 


The mass that flows through the cross-section Q of the tube per second can now 
immediately be determined, namely 


2H 


R 
Q 0 0 





4n n 
(5.19) 


and is proportional to the 4th power of the radius R though the cross-section only 
grows proportinally to its 2nd power. 


Remark: This mass flow was discovered by the prussian civil service hydraulic engi- 
neer Hagen’! by experiments in 1839 without knowledge of (5.18), and independently 





7! Gotthilf Heinrich Ludwig Hagen (1797 - 1884); Berlin 
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also by the french physiologist Poiseuille”? in 1841 when he was interested in blood 
flow. It was also Poiseuille who found formula (5.19). The parabolic velocity distri- 
bution was for the first time mentioned in 1845 by Stokes. 


The experiments of Hagen and of Poiseuille were realized for small tube diameters 
(0.015 - 3.00 mm) and confirmed formula (5.19), thus also the law (5.18) now named 
after them. Reynolds used glas tubes of higher diameters up to 10 mm, he alo filled 
some dyes into the liquids he worked with, and he found out that there is a critical 
value of the parameter which is now called the Reynolds number, and this critical 
Reynolds number is characteristic for the transition of the flow from laminar to tur- 
bulent. He found Reer = 1160, and he recognized that the Hagen-Poiseuille law 
loses its validity in regions of turbulence: The mean flow velocity decreases because 
of a higher flow resistance. 


The drag force of the laminar Hagen-Poiseuille flow in a part of the cylindric pipe of 
length L can be computed by means of (5.14) and of (5.18) in the following way: 
p R 


W =n2TRL (0,u),-p = N2TRL 7 =r R Ap, (5.20) 





where Ap denotes the pressure difference between the two ends of the particular 
part of the pipe under consideration. 


There are also some other cross sections of ducts for which the Navier-Stokes equa- 
tions can be solved analytically under the same assumptions on the flow as in the 
Hagen-Poiseuille case (steady, laminar, incompressible, viscous), e.g. for annuli, 
equilateral triangles etc.”’ 


Situations similar to the flow of a liquid through narrow ducts occur for other flows 
through narrow gaps or slits, e.g. investigation of the flow of lubricating oils in lubri- 
cating research called tribology. 


As already mentioned, there are only very few examples where the Navier-Stokes 
equations can be solved exactly. This led to efforts to simplify these equations wher- 
ever this seemed to be suitable. Some of the important developments concerning 
simplifications in case of small and of great Reynolds numbers will briefly be dis- 
cussed in the next sections. Also for great Reynolds numbers we do not go back to 
the Euler equations but preserve some of the structures of the Navier-Stokes equations 
and study other suitable simplifications, and the no-slip condition will also be kept in 
case of great Reynolds numbers. The case of small Reynolds numbers was particu- 
larly investigated by Stokes and led to a model equation nowadays called the Stokes 


72Jean Louis Marie Poiseuille (1799 - 1869); Paris 
CY. Wang: Ann.Rev.Fluid Mech., 159-177 (1991) 
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approximation (cf. (5.25)), whereas theory as well as practical applications of fluid 
dynamics for small Reynolds numbers was significantly pressed ahead by Prandtl’4 
in his boundary layer theory and his wing theory. 


5.3 Stokes Approximation and Artificial Time 


The difficulties in solving stationary or instationary Navier-Stokes problems or Euler 
problems are mainly due to the nonlinearity of the convection term (u, V) u. 


Another difficulty particularly connected with the case of an incompressible instation- 
ary flow is the fact that the pressure p one has to determine besides the velocity u is 
not represented in the continuity equation by its time derivative, e.g. by an additional 
term «K p (r = const). If this would be the case and if a suitable initial pressure 
po would be prescribed, the instationary problem could be reformulated as a genuine 
initial boundary value problem whose differential equation could be written as 


ape — L[lp|=0 with ẹ = (u, uz, uz, p)” (5.21) 


and could theoretically as well as numerically be treated by means of appropriate 
techniques. 


Let us assume that the solutions of instationary Navier-Stokes problems converge for 
t — oo in each case towards the solution of the correspondent stationary problem 
independent of the initial function uo. 


If this also holds for the instationary problem (5.21) consisting of the Navier-Stokes 
equations and of a continuity equation enriched by a physically unjustified additional 
time derivative of p, the solution of this manipulated problem must coincide with the 
stationary solution of the original Navier-Stokes problem for t — oo because of 


p —>0 for t > œ. 


This leads to the following idea of a method for solving the stationary Navier-Stokes 
problem: 


Introduce an artificial time 7 and enrich the stationary Navier-Stokes equations 
as well as the continuity equation by adding linear combinations of certain terms in 
such a way that after addition of a more or less arbitrary initial function a complete 
instationary initial or initial boundary value problem comes into being whose solution 





1 udwig Prandtl (1875 - 1953); Hannover, Gottingen 
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is expected to have the property to annul the additional artificial terms for T — oo. 
In engineering, this concept is called method of artificial compressibility and is 
occasionally numerically realized. 


Such an artificially enriched stationary problem might look as follows: 


— i i 
ðu (u, V)u+ Vp PE Q, (5.22) 


—0,p = divu 
u=Ħup on OQ, 
u= u", p=po for 7T=0, 


where ull is assumed to coincide with up on OQ. a > 0 and 8 Æ 0 are arbitrary 
constants. They should be chosen in such a way that the numerical treatment of the 
artificially enriched problem is going to become a method as effective as possible. 


Of course, the expected property of the artificially added terms to vanish for increasing 
values of 7 should be guaranteed by a mathematical proof. 


From a numerical point of view, this strategy can lead to an iterative procedure in 
order to solve approximately the stationary boundary value problem in the following 
way: 


Use the more or less arbitrarily added artificial initial function as initial approxima- 
tion, i.e. 


(0] 
Oep et , 5.23 
5 ( A ) (5.23) 


Then discretize (5.22) with respect to the artificial time variable 7 using in every step 
a suitable step size A, 7, e.g. an equidistant step size Ar. 


An example for the semi-discrete approximate equation coming into being is the ex- 
plicit Euler method 


gP =y —a,r-L[p™]  (n=0,1,2,---) (5.24) 


formulated in the style of (5.21). Here, every u!”) is expected to fulfill the boundary 
conditions given in (5.22). yin! will then not only be interpreted as an approximate 
solution of the artificial instationary problem on its nth artificial time level but also 
as an approximate solution of the stationary problem after the nth iteration step. If 
the iteration converges sufficiently fast, it will be stopped after a certain number of ng 
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steps, and glo] (x,y,z) will then be considered to be the final approximate solution 
of the stationary problem. 


At the beginning, this method was used and justified only heuristically by engineers 
who asked the mathematicians for an also mathematical justification. This justifica- 
tion could then at least partially be given by a proof of the fact that in case of the 
Stokes approximation the solution of (5.22) really converges for 7 — oo to the so- 
lution of the stationary problem’>. Here, the truncation errors brought in by additional 
space discretizations are left out of account. 


The Stokes approximation is developed from the Navier-Stokes equations by omission 
of the convective term. This seems to be reasonable in case of incompressible flows 
for small Reynolds numbers. Namely, if (5.11) is written as 


Re- 0; + Re. (ù, V)& + Re-Vp — At = 0, 
one can see that the linear diffusion term dominates the convection term. 


If we omit the tildes and if we restrict ourselves to the stationary problem, after ne- 
glection of the convection term only 


Re-Vp — Au =0 (5.25) 


remains. Finally, we replace Re - p by a new quantity p and enrich the system exem- 
plarily in the way presented in (5.22). The system then turns out to take the form 


du + Vp=AutaViivu } Q, (5.26) 


d-p + B divu=0 





u=0 on ðQ, 
u = ul, p=po for 7=0, 


where also the linear structure of (5.25) was taken into account. Because of this 
structure, it suffices to show the convergence of u for increasing values of 7 to the 
stationary solution in case of vanishing exterior forces and of homogeneous boundary 
conditions. 


The velocity part of the stationary solution of (5.26) then obviously vanishes, i.e. 
u = 0, 


and the following theorem holds: 





St. Fellehner, Dissertation Hamburg 1995. 
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Theorem 5.1: 
The solution u of (5.26)"° converges in the sense of the L2-norm for T — œ to the 
zero function.” 


PROOF: Using the denotations 
Ju-v=((Vu,v),(Vu,v), =, (Vua, »))!, 


(Ju, Jv) = X Oii On, Vi; 
ij 


ae f pgdæ,  lelg=vẸ@,ð in (9), 
Q 

(u,v) = fiu, eae, lül = va in L0), 
Q 


(Ju Jv) = f(u, Jv) da, |J ulo = yJ u, Ju) in L49), 
Q 


taking 
div (J u - v) = (Au, v)+ (J u, J v), 


into account or —by means of Gauss’ divergence theorem— the relation 
Ou 
(Au, v) = ane do — (J u, Jv), (5.27) 
n 
OQ 


and taking finally also into account that the integral along the boundary vanishes, 
(5.26) leads to 


ld T: 1 
SE (h+ illh) = (eu, ue) + slo, Op) 


= (u,-Vp)+(u,Au)+ (u, aV(divu)) 
1 
ee 


= —|J ull — a lldiv ull. 


(p, —G? div u) (5.28) 


7©For convenience, we omit considerations concerning the behaviour of p. 
Here, it is assumed that there is a unique smooth solution of the artificial instationary problem. 
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Thus, 


Hees (iui + zle) 


decreases monotonously, and this the stronger the greater a’® 


lello are bounded so that 


. Hence, ||ul|o and 


Tæ := lim ø(T) 
T—0O 
exists. 
The Poincaré inequality’? then yields in L2“(Q) 


2 2 
lulo < a lJ ullo 


with a constant cı > 0 only depending on the domain Q. Equation (5.27) therefore 
leads to 


d 1 
o'(r) = E (I+ Zll) < -el (5.29) 


where c= 2 > 0. 
Cl 


Formula (5.29) implies 


Habs a lull? d7 (5.30) 


TO 
so that particularly the integral on the righthand side exists, i.e. _||u||o € Lo(R*). 


But also 


(ul) 


is bounded as can be seen by the following arguments: 


Differentiate the differential equation in (5.26) with respect to 7 which leads in a first 
step analogously to (5.29) to 


d 1 
E (loul + ldel) < -elh (5.31) 


*8This can also be used for numerical purposes. 
™cf. e.g. Dautrey-Lions: Mathematical Analysis and Numerical Methods for Science and Technol- 
ogy, vol. 2, p. 126 
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hence to the boundedness of ||3+ul||. Because of 


= (iui) 


2\(u, A-w)| < 2 llull ||O-ullo ; 





the existence of a constant w with 


= (lull) 


follows. The inequalities (5.32) and (5.30) then imply 


<w (5.32) 





jim ||u(7)Ilo = 9. 
Otherwise there would exist an € > 0 and a number 7; > 7 for every 7 > 0 with 
u(r) le > e, and by means of the inequality 


—w (T = T) < lella = lela < w(r — T), 


which follows for every 7 > 7; from (5.32), 
lutra = e- w (rT — r) 


would follow. Thus, also 


Tits 
2 2 
€ We 
a l lullig Ta w 2 w 2w 


would result. 


Now, after the choice of a suitable subsequence of the sequence {J;} the intervals 
(i = 1,2,---) could be assumed to be disjunct and would lead to 


eS aT ta 





SS i I|eu(7) I, dr < œ. 
i m 











This finally would contradict the unboundedness of the partial sums of X` Jj. 
i 
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5.4 Foundations of the Boundary Layer Theory; Flow 
Separation 


Let us now look at the other extreme situation of the Navier-Stokes equations written 
in the dimensionless form (5.11), namely at the case of great Reynolds numbers par- 
ticularly studied by Prandtl. Here, we omit again the tildes used in (5.11), and because 
the Reynolds numbers are great, the behaviour of the solutions can be expected not to 
differ too much from the solutions of the Euler equations®® 


Prandtl published his results in 1904 on the International Congress of Mathematicians 
in Heidelberg®!. Let us now follow his ideas. 


Assume a solid to be put into a viscous flow, let the no-slip condition (1.54) hold 
along the wall IT’ of this solid, i.e. wp = 0, omit the influences of exterior forces, 
and assume the flow to have been irrotational before the solid was put into it. 


Thus, the velocity u could at this time be derived from a potential ® which fulfills 
the potential equation A ® = 0 because of the incompressibility of the fluid. This 
leads to the fact that there was no difference compared to the frictionless situation. 


Prandtl restricted his considerations to a piece I of the surface of the solid which can 
be assumed to be (more or less) plane®’, and this plane was then used as (x, y)-plane 
in an (x,y, z)-coordinate system. In a great distance of the plane, theoretically for 
z — œ, Prandtl expected the flow to keep its character as a potential flow so that the 
Euler equations hold for large values of z. 


The idea that the flow could perhaps be irrotational everywhere outside the solid leads 
for u #0 immediately to a contradiction because the problem 


A® = 0 outside the solid, 


ur = 0 no-slip condition 


has only the solution u = 0. We omit the proof and take into account that a jump 
from the velocity u =O along T` to a potential flow of a velocity u # 0 on the 


8°Tt took many years to make mathematically sure that the Navier-Stokes solutions are close to the 
Euler solutions for small values of 7: K. Nickel showed in Arch.Rat.Mech.Anal. 13, 1-14 (1963), that 
both solutions can under certain conditions and for 7 — 0 asymptotically be looked upon as solutions 
of the Prandtl boundary layer equations (5.48) presented in this section. This also justified subsequently 
Prandtl’s assumption on the small diameter of the boundary layer we are going to use later. 

8! Verh. d. IIL. Int.Math.-Kongr. Heidelberg 1904, S. 484-494, Leipzig: Teubner 1905 

82W, Tollmien showed that all these considerations also hold for curved surfaces provided that the 
curvatures do not vary too much. Cf. Handb. d. Exper.-Physik, IV, part I, 248 ff. (1938). 
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outside of the solid is impossible, too, because the existence of the viscosity term does 
not allow discontinuities of this type. 


Hence, the flow between the surface of the solid and the distant parts of the flow can’t 
be irrotational. The domain of this rotational flow is called the boundary layer. Here, 
it is more or less open what the word distant means in this context from the point of 
view of real applications. In other words: What is the thickness 6 of the boundary 
layer, i.e. from which distance on can the difference between the real flow and the 
original potential flow be neglected ? Of course, the answer depends on a convention. 
Often, 1% difference between the velocities is an accepted value. 


Prandtl assumed the thickness ô to be small compared with the length L used in 
order to write the Navier-Stokes equations in the dimensionless form (5.11), i.e. 


é<1, (5.33) 


and he assumed that the transition from the rotational flow to the irrotational potential 
flow behaves continuously. 


Let all the derivatives occuring in the following considerations exist, let them be con- 
tinuous, and let the dimensionless velocity U of the distant potential flow be parallel 
to T, i.e. parallel to the (x, y)-plane. Hence, 


U =U(z,y,t) = (U(z,y,t) , V(r, y,t) , OF. 


Moreover, the components of U are assumed to be of the same order of magnitude 
as the velocity wo used for the derivation of the dimensionless form (5.11) of the 
Navier-Stokes equations so that U and V are of the order 1 of magnitude. Because 
of the continuous transition from u = (u,v, w)” to U, the assumptions 


u=0(1) and v=O(1) fr z—=ô (5.34) 


with the Landau symbol © are reasonable, and this will certainly also be true along 
all the layers parallel to I, particularly for x,y — 1 where x and y are the 
dimensionless spatial variables in (5.11). 


Concerning the spatial derivatives of the u- and v-components parallel to I as well 
as their time derivatives, Prandtl expected them not to vary very much within the 
boundary layer, i.e. 


Oiu , OiU OW Oyu, Ong s Oyyu , OxV , Oy Ores Oyyv = OU) (5.35) 


If additionally the continuity equation div u = 0 is taken into account, (5.35) yields 
a further result, namely 


0.w=O(1) fr z->06 orfor z,y—>l1, (5.36) 
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so that the no-slip condition leads to 


|w(x,y, z, t)| = fown etae < |O2w ||, Iz, 
0 


w=O(6). (5.37) 
Combining assumption (5.33) with the no-slip condition, the mean value theorem 
leads inside the boundary layer to 

u(x, y, z,t) © O,u(x, y, z, t)z, 
and therefore with (5.34) to 


1 


ozu = O (5) (5.38) 


Analogously also the relations 
1 1 1 
Ov =O (5) 5 zzu =O (=) 5 zV =O (=) (5.39) 


hold. 
We also obtain by the mean value theorem the statement 
w(x, y, 250) = w(0,y, z, t) 0,00 z, y,z,t)-x£x, OR <1; 
and because x can grow up to the order 1, (5.37) yields 
zw = O (ô) (5.40) 
and analogously 


OyW , Ope W , Oyyw = O(S), Aw = O(6). (5.41) 


Moreover, taking (5.36) into account, we additionally find 


ðzzw = O (5) (5.42) 
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Let us now write down separately the three components of (5.11) where we are going 
to note down below each of the terms occuring in the equations their particular orders 
of magnitude as far as they are already known. Here we follow a representation given 
by H. Schlichting in the first edition of his book on boundary-layer theory:83 


Ou + u Ozu + v Oyu + w Ozu -k (rcu + Oyyu + zzu) 


1 1 1 1 1 5 + 1 1 2 
Ov + u sv + v Oyv + w av -Rz (brzu + Oyyv + zzv) 
1 1 1 1 1 ô 3 1 1 z 
Ow + u rw + v Oyw + w ðw -K (rew + Oyyw + zzw) 

6 1 ô 1 ô ô 1 ô ô + 
(5.43) 


Of course, the pressure gradient parallel to the plane is assumed to be bounded in all 
real applications, i.e. 


ozp ) oyp = O(1) : 


But this implies that the first two equations can only be valid with respect to the order 
of magnitude if 


1 E 1 
Re = O (=) , ie 6 =O (=) : (5.44) 


and this leads because of the third equation to 





zp = O(ô). (5.45) 
The equation (5.44) shows that Prandtl’s assumption on a small thickness of the 
boundary layer for great Reynolds numbers does not lead to a contradiction. 


It also shows that the boundary layer theory can only be accepted as an approximation 
for the Navier-Stokes equations in the case of great Reynolds numbers. 





H, Schlichting: Grenzschichttheorie; Karlsruhe, Braun 1951 
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Because of 
zZ 
|p(z, Y, Z, t) a plz, Y, 0, t)| = | O-p(x, Y, Ç, t) dÇ < |z| Il zp lloo ’ 
0 

the pressure differences in z-direction within the boundary layer are of order 6°, 
hence very small for great Reynolds numbers. We therefore put 

Ozp =0 ’ 
so that p does approximately not depend on z within the boundary layer: 


p= p(z,y,t). (5.46) 


The pressure gradient in the boundary layer thus coincides with the well known pres- 
sure gradient of the potential flow outside the boundary layer where the Euler equa- 
tions hold. Using dimensionless variables and taking W = 0 into account, we 
therefore obtain approximately 





OV +U0;,V+V0,V = —Oyp 
so that O,p and Oyp are known quantities. 


If now all the terms in (5.43) of order 6 will be neglected compared with the terms 
of order 1, the system of three equations reduces to only two equations, and these 
equations read in dimensionless form as 


Ou + uru + vdyu + wd,u — vdz,u = ——Ozgp 








Ov + udgv + vdyv + wd,v — vdz,,0 = ——Oyp 


(5.48) are the Prandtl Boundary Layer Equations added by the incompressibility 
condition 


Ozu + Oyv + O,w = 0. (5.49) 
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The velocity component w still remains in the system (5.48) but does now only play 
the role of a parameter which will later be fixed by (5.49) and by boundary conditions. 
An example of such boundary conditions reads as 








u v w 0 for z=0 
u U(x,y,t) 
v |=| Vig) for z>. 
w 0 


z — œ and ô < L seem to be contradictory but one has to take into account that the 
boundary layer is of an infinite thickness from a theoretical point of view, and only 
their effective thickness judged from a physical point of view is very small. 


If the boundary layer is affected by sucking off or blowing up parts of it, the boundary 
conditions have to model these facts, e.g. by replacing the equation w = 0 in (5.50) 
by 


w = wolz,y,t) for z=0, 


provided that the direction of the sucking or blowing coincides with the w-direction 
and that its flow velocity wo is of order O(ô), i.e. 


a-o (ga) 


Otherwise some of the terms expected by Prandtl to be very small can no longer be 
neglected. 





Besides boundary conditions also initial conditions have to be formulated, and this 
not only with respect to time in case of instationary problems but also in every case 
with respect to space because the Boundary Layer Equations are of parabolic char- 
acter. This parabolicity often includes simplifications as far as the construction of 
approximate analytic solutions of the Navier-Stokes equations is desired. But also 
with respect to numerical purposes Prandtl’s theory can often make life easier. 


If the problem to be treated concerns an instationary flow, the task to begin with is the 
formulation of initial conditions with respect to time, i.e. conditions of the type 


u(x, y, z,0) = fila, y, z) 
v(z, y, z,0) = fo(x,y, z) ’ 
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and in every case also with respect to space, i.e. the stipulation of suitable initial 
profiles of the u- and v-components of the flow velocity. 


There is no general consensus how the initial profiles should look like, but it should 
be noticed that the influence of the initial profiles decreases for increasing distance 
from the area where these profiles were prescribed. This follows from the parabolic 
situation. 


These spatial initial conditions may take into consideration that a boundary layer 
comes into being only after a certain distance from the edges of the solid put into 
the fluid, e.g. from the front edge of a wing: 


u(Xo, Y, z, t) = gly, z, t) 
v(xo, yY, 2, t) = gi2(y, Z,t) 


and 


u(x, Yo, z,t) = gaia, z,t) 
(a, Yo; z,t) = g22(x, Bt) : 


Of course, the conditions prescribed should go well together, e.g. 


fi(xo,y, z) = gıı (y, 2,0) (5.54) 


etc. 


Many phenomenons of flows taking place along the surface of a solid can quanti- 
tatively be explained by means of the boundary layer concept. One of them is the 
separation of the boundary layer flow from a curved surface. 


Here, we are going to describe this separation briefly and only qualitatively in case 
of stationary 2-D flows: Assume an ideal fluid to flow along a curved surface. A 
lateral flow with main flow velocity u = (uə,0, oy! , Ugo > 0, around a circular 
cylinder of infinite length whose axis coincides with the z-axis is a 2-D example. The 
velocity of particles arriving at the surface of the cylinder increases due to Bernoulli’s 
equation (1.31) because they move at first against decreasing pressure. The velocity 
then begins to decrease because the pressure grows again. Finally the particles move 
behind the cylinder with their original speed (cf. Fig. 26). 


In the case of a viscid flow, a certain amount of the kinetic energy of the particles will 
be annihilated in the first phase by friction so that there is in the second phase not 
enough energy available to bring the particles back to their former speed. It can even 
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Figure 26: Ideal flow around a cylinder 


happen that the particles stop and that the velocity reverses its direction because of the 
increasing pressure, i.e. separation of the boundary layer takes place. 


In the 2-D situation, the point of separation at the curved surface which is in two 
dimensions a curve I’ can be well described mathematically. For this purpose, let u 











a 


Figure 27: Flow with separation 
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be the particular component of the flow velocity parallel to T and let ôu be the 
normal derivative of this component. The separation point P € I is then obviously 
given as the point where this normal derivative vanishes (cf. Fig 28): 


nul P) =0. (5.55) 





T 


Figure 28: Explanation of (5.55) in case of a 2-D boundary layer flow 





Figure 29: Boundary layer separation on a wing 


Beyond the point of separation, the flow loses its buoyancy effect. Hence, one tries to 
design a wing in such a way that the point of separation is situated as near as possible 
at the back edge of the wing, e.g. by sucking off the boundary layer. 


It should be mentioned that there is no general mathematical definition of separation 
in case of 3-D flows. 
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5.5 Stability of Laminar Flows 


The beginning of the transition of a flow from the laminar state to a turbulent one can 
be looked upon as a transition from a stable behaviour of the solution of the Navier- 
Stokes equations belonging to this flow to instability. 


In order to approach this idea, let us for convenience consider a situation easy to 
survey, namely an incompressible 2-D flow and the knowledge of a stationary average 
laminar flow in z-direction whose velocity U only depends on the second space 
variable y: 


The pressure of this flow is denoted by P(x,y), and (U, P) is assumed to fulfill the 
given incompressible Navier-Stokes problem** and to be one part of a superposition 
whose other part consists of infinitesimal incompressible disturbances (ù, p) with 


oa) ETA 


This second part is not necessarily a statinary flow and, besides (U, P), also the 
complete flow 


u=U+u, p=P+p (5.56) 
is assumed to be a solution of the given Navier-Stokes problem®?. 


Because the disturbances are infinitesimal, quadratic terms can be neglected compared 
with linear terms, so that putting (5.56) into the Navier-Stokes equations (5.10) and 
into the continuity equation then yields the relations 





2. 
Depo EE He PP GES (SS +i) 
dy p y 














p dy? 
~ x 1 ch agen 3 
OO + U zð —OyP + —Oyp =vAd (5.57) 
p p 
rù + Oyo = 0. 


8t is absolutely problematic to suppose such an average flow to exist. 
This assumption is called the closing, and the treatment of the instability problem in this way is 
called the Reynolds averaging. 
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Because the average flow solves the Navier-Stokes equations on its own, (5.57) re- 
duces to the linear equations 
dU 


1 
O+ Uù +¢60— 4+-0,p=VAt 
dy p 


1 
HO+U Azo +—Odyp =v Ab (5.58) 
p 


A, t + yo = 0. 


If the first equation in (5.58) is differentiated with respect to y and the second one 
with respect to x, the pressure can be eliminated so that there remains only one 
equation for u and v besides the continuity equation. The no-slip condition along 
surfaces of solids leads to additional boundary conditions. 


Let us now suppose that the disturbances result from the superposition of waves propa- 
gating in x-direction where each of them fulfills (5.58). These waves can be expressed 
by 


< p(y) ) i(ax—4t) 
u= e 5.59 

( vy) ae 
and are called Tollmien-Schlichting waves. Here, a € C, ~y € C. The particular 
case a € R describes the investigation of the temporal stability, whereas an inves- 
tigation of the case a € C, ~y € R concerns the so-called spatial stability not 
discussed here. 


The continuity equation leads to 


and the equation that remains after the eliminition of p then yields 


(U — a)" ap) — UY = -E (YO — 202" + ay). (5.60) 


The linear homogeneous 4" order ordinary differential equation (5.60) for ~ with 
complex coefficients is called the Orr-Sommerfeld equation. Boundary conditions 
are for instance the no-slip condition 7 = 0 for y = 0 and the disappearance of the 
disturbances for y — oo. 


Thus, for every fixed a € R, an eigenvalue problem arises where the complex phase 
speed of perturbation y plays the role of the eigenvalue. The parameter v*® is 


; ; at 7 1 
86 or, if the problem is formulated in dimensionless form, Re 
e 
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considered to be already known from the laminar average flow. This laminar flow 
corresponds to the trivial solution ~ = 0 of problem (5.60), and one asks for 
nontrivial solutions bifurcating from it. 


If the sign of the imaginary part of y is negative, the disturbance decreases when 
time grows, and its amplitude does not depend on time if this imaginary part vanishes. 
But an unstable situation occurs for 


sgn (Im 7) > 0. (5.61) 


The Orr-Sommerfeld equation can under realistic boundary conditions for U nor- 
mally not be solved explicitely though it is linear and homogeneous. Numerical pro- 
cedures known from numerical bifurcation theory have to be used. 


The smallest value of the Reynolds number giving rise to the validity of (5.61) is 
called the critical Reynolds number Re, and is of interest because it is the smallest 
Reynolds number for which the flow can become unstable. The curve sgn (Im y) = 0 
obviously separates the area of instability in the (a, Re)-plane from the area of a stable 
laminar flow and is called the curve of neutral stability. 


It should be noted that an instability result due to this linearized theory of disturbances 
does not necessarily ensure that the real underlying flow can already be considered to 
be completely turbulent. A measure accepted by engineers is the so-called e?-rule: 
The flow is really a turbulent one if the small amplitude of the initial disturbance has 
increased by the factor e? ~ 8100. 


6 Existence Proof for Entropy Solutions by Means 
of Discretization Procedures 


Though we have already studied some properties of weak solutions of conservation 
laws, we did not yet proof the existence of this type of solutions, particularly of en- 
tropy solutions. The only exceptions were some scalar problems whose solutions 
could explicitely be specified. 


We are now going to attack the existence problem. 


6.1 Some Historical Remarks 


Existence proofs for solutions of ordinary or partial differential equations are often 
structured in the following way: In a first step, the differential equation will be dis- 
cretized as normally done in Numerical Analysis in order to create numerical proce- 
dures. Thus, a finite dimensional approximate version of the original problem solvable 
on a certain grid for every suitable step size parameter h comes into being. If the 
solution of the discrete problem is only given on a set of discrete points, it will be ex- 
tended to the inter grid points by interpolation or other techniques in such a way that 
it belongs to the particular function space the solution or weak solution of the original 
problem is expected to belong to. This extension is called reconstruction (cf. the 
special case (4.2)). Moreover, we assume that all the computations to be made were 
theoretically realized without occurence of round off errors. 


If it can be shown in a second step that a sequence £ of these reconstructed approx- 
imate solutions determined for every step size parameter hj of a sequence {h;} of 


positive step size parameters with lim = 0 is compact in the function space under 
jroo 


consideration, convergent subsequences of £ exist. It can often be shown in a third 
step that under certain conditions the limits of those subsequences solve the original 
problem. Thus, solutions exist in this case, but different convergent subsequences can 
lead to different solutions of the original problem because uniqueness was not yet 
ensured. 


Mathematical Models of Fluiddynamics. Rainer Ansorge 
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A well known example of an existence proof of this type is the particular proof of 
Peano’s theorem on the existence of a solution of the initial value problem 


y' = f(x,y), y(zo) = yo, f continuous 


where Euler’s polygon method is used in order to construct approximate solutions 
in the space of functions continuous on a compact interval and equipped with the 
maximum norm. 


Another example is the existence proof of Courant, Friedrichs, Lewy where this type 
of proof was for the first time applied on initial value problems of systems of quasilin- 
ear hyperbolic partial differential equations. The proof was published in their famous 
paper of 192837, and it was in this paper that a condition concerning the step ratio was 
formulated in order to ensure numerical stability of the numerical procedure, later 
called the CFL condition. This condition ensures that the region affecting the solution 
at a point Q where the solution is expected to exist is covered by the region that af- 
fects the approximate solution at Q, and this condition is extremely important when 
treating conservation laws by means of explicit numerical methods’. 


Of course, because Courant, Friedrichs and Lewy asked for smooth solutions, they 
could only give conditions for local existence. 


Under certain restrictions concerning the initial functions, J. Glimm gave an existence 
theorem for global weak solutions of conservation laws in a paper that also became 
famous®?, and one of his tools was also a certain finite difference procedure. 


We are now going to study a general theory on discretizations of partial differential 
equations that can also help to prove existence theorems for conservation laws as 


occuring in fluid dynamics”. 


6.2 Reduction to Properties of Operator Sequences 


In 1948, Kantorowitch published his important paper on the connections between 
91 


numerical mathematics and functional analysis”. 
Particularly in case of linear partial evolution equations, Lax and Richtmyer provided 
in 1956 a functionalanalytic frame for the construction and investigation of finite dif- 


87Math. Ann. 100, 32-74 (1928) 

88For more details cf. Section 7.3. 

8° Comm. Pure Appl. Math. 18, 697-715 (1965) 
R, Ansorge, ZAMM 73, 239-253 (1993) 
*!Uspeki Mat. 3, 89-185 (1948) 
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ference methods” (cf. section 7.4). Their theory was later generalized with respect 
to different aspects. 


Inspired by the Lax-Richtmyer Theory but besides the particular question of the dis- 
cretization of differential equations, Stummel published in 1972 his general theory 
on discretization algorithms”*. Though it was already rather general, it was still influ- 
enced by linear structures underlying the given general problem. This theory was later 
also generalized into different directions and applied to several types of problems such 
as ordinary or partial differential equations, integral equations, optimization problems 
etc. Since then, this theory was called Functional Analysis of Discretization Algo- 
rithms (Vainikko). 


But not yet included into this theory were weak solutions of differential equations or 
problems whose solutions are not necessarily unique. 


Let us therefore in a first step try to describe an imbedding of discretization techniques 
for problems of this type into the general theory and look for problems given in the 
following general form: 


Let X, Y , Xn , Yn be topological, usually metric spaces, and let us ask for an 
element v € X which fulfills 


Av=a, A: X ~Y a given operator. (6.1) 


Let the numerical procedure constructed in order to solve this problem also be de- 
scribed in a general form by 


AnUn = ün An : Xn —7Yn, nEN, (6.2) 


where {an|an € Yn, n € IN} is a given sequence approximating a € Y ina sense 
still to be specified. 


Neither Xn C Xnii nor Xn C X and neither Yp, C Yny1 nor Yp CY (n= 
1,2,---) are assumed to hold. 


Assume that there are either mappings 
Pn: X > Xn, dn: Y — Yn 
Pn | qn | 
Xn An Yn 


°2Comm. Pure Appl. Math. 9, 267-293 (1956) 
Proc. Conf. Num. Anal. Dublin 1972 (ed.: J. Miller), 285-310, Academic Press 1973 
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or an imbedding of the numerical procedure, just for theoretical purposes, into the 
frame of the original problem. This can for instance be done in the sense of 


Xn Cx, Yn =Y (n=1,2,---), 
but not necessarily 


Ree ay Le 


Just for convenience, let us suppose such an imbedding to be realized, but the given 
original problem does not necessarily need to be linear, and its solutions are allowed 
to be weak in a sense to be specified. 


Let the original problem be described by 
Av=w, (6.3) 
with A: X >Y, X CX. 


The weak formulation of (6.3) which then will be the problem really to be solved is 
supposed to be described in the following way: 


Let J be an index set; let {a(®),@ € J} C Y bea given set and {A(®),® € J} 
a set of operators with joint domain D CX, D-Y. 


We ask for an element v € X with 


A(®)v = a(®) y Oey: (6.4) 


Definition: 
The elements of 


O := {v € X|v solves (6.4)} 
are called weak solutions of (6.3) if the following implications hold: 
If v solves (6.3), then v € © ; if v € X NO, it also solves (6.3). 
The elements ® € J are called test elements. 


Let the numerical method be described in an also rather general form by 


AnUn = Gn (n = 1,2,---) (6.5) 
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where Ân: X, > Z (n = 1,2,---) with a certain topological space Z. Here, 
{ân} C Z is a given sequence compact in Z. 


Assume 
On y= {un E€ Xn|Un solves (6.5)} £ 0 (n = 1, 2, _ =) , 


but every ©, may contain more than one element. The solutions v,, of (6.5) are 
called approximate solutions of (6.4). 


The numerical method (6.5) is assumed not to depend on test elements because com- 
puters cannot understand what test elements are. But we suppose that also this method 
can be formulated weakly by means of a sequence of operators 


{An(®)|®@ € J, nE IN} with A,(®): Xn > Yn, 
and of a sequence {an(®)} with 


lim an(®) = a(®) for every fixed ® € J, 


n—> oo 


so that 
An(®)un = an(®) VPEJ, Vin € On (6.6) 
holds. 
(6.6) is called a weak formulation of the numerical method. 
Some concepts help to prove convergence results: 


Definition: 
A pair |{Cn} , C], consisting of an operator sequence {Cn} and of an operator C is 
called asymptotically closed if the implication 


Un >v A Chtpnrz > Cv=z 


holds. 


Definition: 
An operator sequence {Cn} is called asymptotically regular if the implication 


{Cnn} compactinY = {un} compact in X 


holds. 
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Definition: 
Method (6.5) is called convergent if set convergence 
0, —-O 
is ensured in the following sense: 


{On} is discretely compact, i.e. every sequence {Vn|Un E On ;n =1,2,---} is 
compact in X, and if v is a limit of a convergent subsequence, v € © follows. 


Remark: 
O, — © implies © ¢ Ø. In other words: There is a solution of (6.4), i.e. a weak 
solution of (6.3). 


6.3 Convergence Theorems 


We can now state a convergence theorem: 
Theorem 6.1: 
(i) Let [{An(®)} , A(®)] be asymptotically closed for every fixed ® € J and let 
(ii) {An} be asymptotically regular. 
Then: 


0, — O. 


PROOF: {ân} compactin Z => { Anvnlen € en} compact in Z, and this for 


every sequence {Vn|Un E€ On}. 


Thus, each of the sequences {vn} is compact in X because of assumption (ii), i.e. 








J {unun € On ,n € IN’ CN} , JveX with uy ov (neN'’). 
Here, IN’ and v are independent of ©. 
But because vp, fulfills also An(®)vn = an(®), and because of 


anl)  a(b) VOET, 











supposition (i) yields A(®)v = a(®), i.e. v € ©. 





4A proof of ©,, — © is therefore also an existence proof. 
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If one is only interested in a particular problem, i.e. in a particular right hand side of 
equation (6.3), assumption (ii) can obviously be replaced by the weaker requirement 


(ii*) {On} discretely compact. 


If uniqueness of a weak solution is not guaranteed, sometimes additional restrictions 
are formulated in order to pick out of the set of weak solutions the particular one 
relevant from the user’s point of view. And as in the case of scalar conservation 
laws, this restriction often takes the form of an inequality the relevant solution has 
additionally to fulfill. Therefore we are going to call such a restriction generally an 
entropy condition, and we formulate this condition in a general setting as 


B®)v<0 VOeES, (6.7) 
where {Bô EJ \ denotes a set of nonlinear continuous functionals mapping 
X into R, and where J is a certain index set which may differ from J. 

We suppose that there is a uniqueness theorem available”: 


There is at most one solution of (6.4) that fulfills also (6.7), called the entropy 
solution. 


If the original problem is completed by an entropy condition in order to characterize 
the relevant solution, it seems to be reasonable to discretize also this condition by 


B (jun <0 VGES (n=1,2,---), (6.8) 


and to accept only such approximate solutions vn € S» which additionally fulfill 


this discrete entropy version. Herewith, the functionals B,(®) are also assumed to 
be continuous. 


We suppose 
{(6.5), (6.8)} has at least one solution ên € On for every fixed n € N. 


(6.8) is called a sequence of discrete entropy conditions, but we do not expect a 
discrete entropy solution 


n (n € N arbitrary) 


necessarily to be unique. 





Theorem 3.3 is an example. 
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Definition: 
An operator sequence {Cn} is called continuously convergent to an operator C if 
the implication 


Un > V => Chin > Cv 
holds”©°’, We then write 

Cee: 
We can now supplement Theorem 6.1 with another important remark: 
Theorem 6.2: Let 

Ôn := {in € Onlôn solves (6.8)} Æ 0, 
and let 

Bn (Ê) => B(®) be fulfilled for every fixed ® € J. 
Moreover, assume the suppositions of Theorem 6.1 to be fulfilled. 
Then there is an entropy solution © of {(6.4),(6.7)}, and 

Ô, — {ô}. 
PROOF: 


Because of the validity of Theorem 6.1, i.e. because of O,, — O, there is a subse- 
quence 


{in| in € On, ne W cm} 


and a weak solution v € © with {ôn} >v, ne WN’. 


The continuous convergence therefore yields By(®)in > B(®)v. 


From (6.8), B(®)v < 0 follows. Consequently, v is an entropy solution, and 
because the entropy solution was assumed to be unique, v coincides with this entropy 
solution, i.e. v = 8, and not only a subsequence converges to 0 but the full sequence 
converges: 











lim n = Ê. 





%P, du Bois-Reymond: Sitz.-Berichte, Preussische Akad. d. Wiss., 359-360 (1886) 
"R. Courant: Göttinger Nachr., 101-109 (1914) 
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Conditions for the property of continuous convergence on metric spaces follow from 


a theorem due to W. Rinow”? : 


Theorem 6.3: 

If the operators C , Cn (n = 1,2,--+) mapping a metric space X into a metric 
space Y are individually continuous, then Cn —> C' is equivalent to the simultane- 
ous validity of the two statements: 


(j) {Cn} is a sequence of equicontinuous operators , 


Gj) Cn > C holds pointwise on a subset X eee 


6.4 Example 


As an example, we treat the scalar version of problem (2.1). Thus, let x € R be 
a spatial variable and t > 0 a time variable, and let Q be again the upper (x, t)- 
halfplane. 


We ask for solutions v : Q — R of the scalar conservation law 


w + sfv) =0, t>0 
(6.9) 
v(x,0) = v(x) with vp € BV(R)N L(R), 


where the flux f € C! is strictly convex and has its minimum at zero which can be 
assumed to be zero: f(0) = 0. Let 


I'o = max {|f"(v)], |v] < llvolle } <% - (6.10) 
Because of possible discontinuities we are interested in weak global solutions v € 
X := L of 

A(®)u := -f [ (v2 toaa) dx dt 

Q 
(6.11) 


= f 22.0) dx = 0 VO ey = CEN). 
Let 
S=S(v,c), Rx ROR 


*8Die innere Geometrie der metrischen Räume, p. 78 ff., Berlin-Géttingen-Heidelberg: Springer 1961 
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be a one parameter set of so-called entropy functions continuously differentiable for 
every fixed c € R with respect to v, possibly with an exception at v = c. 
F =F(v,c) (Rx R-— R) is called the entropy flux, if 

O,5(v(a,t),c) + OF (u(x,t),c)=0 on Q 


holds for every fixed c € R and for every pair (v , Q’), for which v is a weak 
solution on Q but even a smooth C1-solution on X C Q. 


If a solution v, which is weak on Q, fulfills on Q additionally the inequality 
O,5(v(2,t),c) + Op F (v(z,t),c) <0 Vee R 
in a weak sense, i.e. 


- | | {ae@.n8o(e,,<) + dle, t) Polet), o} da di 
= f Eei) dx <0 6.12) 


VYceER , V nonnegative P E€ J, 
then v is called the entropy solution, denoted by 0. 
Thus, test elements in this example are the pairs 
Ê := (,c) with Ê c€ J:= {(8, c) EJ, 8 >0,ceER}, 
and the integral inequality can therefore be written as 
B(w <0 VOeT (6.13) 
in accordance with the general theory. Particular choices of ĝ in the literature are 


ï= S(v) , S independent of c but strictly convex, 
b={GeJO>0}cJ 
(Lax; cf. (3.11)) , 
S(v,c) = |v — cl (Kruzkov; cf. (3.44)) , 


p = > 
S(v,c) = { a 0 E a Me (Harten, Hyman, Lax”) 


Comm. Pure Appl. Math. 29, 297-322 (1976) 
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Each of these choices then implies a particular choice of F according to (3.13). 


Let the numerical procedure used in the general setting be put in concrete form by a 
(2k+1)-point finite difference scheme: 


v+1 v — qv 
a aot a =0, (6.14) 
where v counts the number of time steps, j the number of spatial steps, starting from 
j+1/2Az 
v? = — Uo(x) dz. 
(j-1/2)Ax 


Here, the numerical flux 
— V V V V V 
Ij42 = 9(Vj_ Rp Us_k4a0° °° Vj Uj449°°° Utk) 
is assumed to be Lipschitz continuous and consistent with the original problem, i.e. 


glw,w,:--,w) = flw), YWER. (6.15) 


Methods of this form are called methods in conservation form. 


We restrict the choices of the step sizes by the requirement that the CFL condition, al- 
ready introduced in the beginning of this chapter, has to be respected. For the problem 
under consideration this means 


At 
0 < à := — = const < 


1 
Dz FR (6.16) 


Let Ax = O(4) (n € N), and put 


(=e @ <Gta)Ae 


Un(z,t) = vj for { yAt< t <(v+l)At 





(G= Oe ee) ee Ti 2 











If Xn is then chosen as the space of functions defined on Q and constant on these 
rectangles, the finite difference scheme can be formulated as 


Antn=0 (n=1,2,:::) (6.17) 


with operators An : Xn 7 Z:=R. 
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Multiplication of (6.14) by test functions ® € J and integration leads to 


An(®) un = 


(v-+1) At (+5) Aa) 


2 2 J | RN jet — Uy +A (t, = g) dxdt = 0 
oa v At (j-4) Az) 


and describes the operators A,(®) as well as the right sides a,,(®), which equal 
zero in this context. 


Using the abbreviation 


Sete) — S¥(e) ; Gle) — G_a(¢) 


< a 
x x <0 VceR, (6.19) 


where 
v _ v yov v 
Gja (c c) = Gluj_egi tt 9 U5) U5410° °° ,Uj+k C) 


denotes a Lipschitz continuous numerical entropy flux to be specified and consistent 
with the entropy flux F in the sense of 


G(w,w,::- ,w,c) = F(w,c) V(u,c) € R?. (6.20) 


For k = 1 (3-point case) the flux splitting choice 
G(a, b,c) = F+ (a,c) + F_(B,c) (6.21) 
with 


O B20 + ROOT Bag | Bx0 


is one of several possible choices of G, where F is known as soon as S'is chosen. 
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The weak formulation of the discrete form of the entropy condition needed for the 
use of Theorem 6.2 arises from multiplication of (6.19) by nonnegative test functions 
® e J and integration: 


Bos: = [ fee + [S(un(a,t + At), €) Sonta, t), 0) 
+ = E A ee PA (6.22) 


— G(vn(a — Az, t), un(z,t), c)] dx dt <0 


V@eJ,ie 
t) — At) ~ 
= JE (2, Hee re) S(un(a,t), c) dz dt 
At —oo 
At aj+Ax/2 


</> f P(x, t) da S(v? ,c) dt (6.23) 


0 . a,—Ax/2 
O(x + 42,t)- O(a - 22,4 
-f [EH tonto 48.0, vn(xt+A2,t), 0) dxdt 
Q Ax 


<0 YyYĝcĪ. 


Thus, B and B,, from (6.7) and (6.8), respectively, are put into concrete form provided 
that a 3-point scheme is under consideration. If then, for instance, the Harten-Hyman- 
Lax choice of Š is taken into account and G is chosen according to (6.21), then 
B(®) in (6.12), (6.13) is continuous for every fixed ®. The functionals B,(®) (n = 
1,2,-++) in (6.22) are equicontinuous, and pointwise convergence B,(®) > B(®) 
can easily be shown provided that the joint domain of these functionals is restricted to 
elements v € X* ; herewith, X* is the space of functions v € X restricted to the 
domain supp ® which is equipped with the L1-metric on X*. 


Hence, the assumptions (j) and (jj) of Rinow’s theorem are fulfilled, i.e. 
B,(@) — B®), 


and Theorem 6.2 applies as far as Ôn egy 


100 





cf. footnote 90 
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The explicit!°! Engquist-Osher 3-point flux splitting scheme!, characterized by 
(6.14) with 


Ged = guj, v1) = f- (vy) + £4 (7) (6.24) 


2 


0 for a<0 (a) for a<0 ’ 


pea for a>0 . E. fr a>0 


together with the Harten-Hyman-Lax choice of S, is an example. (6.24) then de- 
scribes the operators An (n = 1,2,---) in (6.17)!%, hence also the operators 
A,(®) via (6.18). The method simulates an important property of weak solutions 
of scalar problems of type (2.1) provided the CFL condition is fulfilled, namely their 
monotonicity, hence its Total Variation Diminishing or TVD property and, there- 


fore, its Total Variation Boundedness or TVB. 
We interrupt the proof of our convergence theorem for some explanations: 


A method used for the numerical treatment of scalar conservation laws of type (2.1), 
represented in the general form 


n+1 _ n n n 
Ug g= H (v_p » UJ_k419 a Ui+k) ) 


is called monotone if the partial derivatives of the function H with respect to all its 
variables are nonnegative. 

This imitates the following property of the weak solutions of (2.1): If two initial func- 
tions ŭo and vo fulfill the inequality 


vo(x) > volz) YTTER, 


this inequality holds everywhere in Q for the weak solutions (x,t) and v(x,t) 
belonging to these initial functions, respectively. 


TV-boundedness in this context (and formulated not only for scalar equations but 
also for systems) means that X`; \|\Vj.. — V4|| K for all n, and TVD meth- 
ods introducded by Harten'™ simulate the fact that the total variation of weak solu- 
tions of scalar problems does not increase, and it is obvious that TVD implies TVB. 
The fact that monotonicity of a scheme implies its TVD property was shown in an 


Rm Ś 


''Tt follows from this explicitness thatO, 40 , (n=1,2,---). 
'©Math.Comput. 34, 45-75 (1980) 

3a,=0 , (n=1,2,---). 

104STAM J. Numer.Anal. 21, 1-23 (1984) 
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appendix to the paper of Harten, Hyman, Lax (cf. footnote 98), written by Barbara 
Keyfitz. 


We continue proving the convergence theorem: 
It was LeVeque!© who showed that ©,, is discretely compact! provided that vo is 


of bounded variation. It is here where we put our requirement vp E€ BV(R) to use. 
Thus, assumption (ii) of Theorem 6.1 is fulfilled because (ii*) is fulfilled. 
The property (i), i.e. 

[{An(®)}, A(®)] discretely closed , 


remains to be shown in order to show the convergence of the method by means of 
Theorem 6.1. But this can in case of TVB methods for scalar equations immediately 
be deduced from the Lax-Wendroff Theorem!” . 


This ends the proof for the existence of an entropy solution of our scalar conservation 
law example. 


As a particular result, we notice that a convergent monotone consistent scheme in 
conservation form converges to the entropy solution of the original problem. 


The Lax-Wendroff Theorem refered to in the proof of the convergence theorem reads 
as follows: 


Theorem 6.4: 


Let X = at = const and let At = O (+) (n € N) for n > œ. Define 
V” € L to be the step function 


Ax Ar 
V"(a,t) = V} for ty oy Sesika , tye St Z tras 














(j = Oya aed ee ; n=0,1,2,---) on Q= {(x,t)|xrEeR,t>0}, 
where the numerical values Ve were computed by means of an explicit TV-bounding 


finite difference one-step scheme whose numerical flux g is consistent and Lipschitz 
continuous and which starts from 














aj+42 
1 . 
Ve / Vo(€)d€ (j =0,£1,42,---). 
Ax 
Tj- 


105Numerical Methods for Conservation Laws, p. 164. Basel-Boston-Berlin: Birkhäuser 1990 
106 though he used another terminology 
107 Comm. Pure Appl. Math. 13, 217-237 (1960) 
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Assume that V” converges for n — œ, i.e for At — 0), to a function V in the 
L'°°-topology, i.e. with respect to the Lj-norm on each compact subset of Q. and for 
each of the components of V”. 


Then V is a weak solution of (2.1). 


We omit the proof of the Lax-Wendroff Theorem but notice that this is not yet a 
convergence theorem! The theorem only says that the limits of convergent sequences 
of approximate solutions —as far as they do exist— are weak solutions of the original 
problem. 


In order to finish off a proof showing the important property of a method to converge, 
one has to look for properties which ensure the validity of the assumptions of Theo- 
rems 6.1 and 6.2 (as it was realized in the case of the Engqist-Osher scheme). 


It should be mentioned that the general convergence theorem of this chapter does of 
course not only apply in case of initial value problems for scalar conservation laws but 
also in case of certain nonlinear elliptic boundary value problems numerically treated 
by means of projection methods etc.!° 





108R, Ansorge and J. Lei, ZAMM 71, 207-221 (1991) 
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7.1 Some General Remarks 


We want to describe the principal possibilities of discretizing hyperbolic conservation 
laws and describe their common features as well as their differences!’. For the sake 
of simplicity we consider the generic scalar equation 


L(u) := bju + O,u+ Oyu = 0. (7.1) 


There are three basic classes of discretization methods for (7.1), namely the Finite 
Element Method (FEM), the Finite Volume Method (FVM), and the Finite Differ- 
ence Method (FDM) already introduced briefly in section 4.1 and already used for 
more theoretical purposes in section 6.4. Although these names sound quite definitive 
there are very many particular schemes behind each of the three classes of methods. 


In the FDM we begin by discretizing time 
G” := {nAt | At > 0,n € NU {O}} 
as well as space 
Gij := {(iAz, jAy) | Av, Ay > 0, i,j € Z} 


and replace the differential operators in (7.1) by difference operators, for example: 


Ax, jA DAŁ) — uliAr, jAy,nAt 
dite phy Ah S MAN As (nF DAD An Jun) 
wj UR 

At 





10° cf. also Chapter 6 


Mathematical Models of Fluiddynamics. Rainer Ansorge 
Copyright © 2003 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim 
ISBN: 3-527-40397-3 
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and 


Az, jAy, nAt) — i — DAZ, jAy. nAt 
O0,u(iAa, jAy, nAt) x uliAx, jAy, nAt) ~ ulli = LAs, jAy, nAt) 


Aa 
= Uj ina 
AT 
Oyu(iAa, jAy,nAt) x WR An U Den nat) 
= Us — Uj- 
Ay i 


Replacing the derivatives in (7.1) by the differences leads to the difference equation 


n+1 n n n n n 
iiaa 2 sg dg a y 
At Aa Ay , 


where the superscript h indicates a discrete value (very often one encounters Ax = 
Ay and calls this value simply h). Since we decided to approximate the spatial deriva- 
tives at time level nAt the resulting difference equation is explicit, i.e. the solution at 
the resulting time level (n + 1)At is explicitly computable through 


Uij SUL; Ae (up, = ugay) T Ay (uf, = Uy 5-1 


If we had decided to take the spatial differences at t = (n+1)At we would have ended 
up with an implicit equation, requiring the solution of a system of linear equations 
at each time step. Note that in (7.2) we still have used the exact solution u of (7.1). 
Since there will be errors due to the replacement of differential quotients by difference 
quotients we can not assume (7.2) to hold for u. Thus, (7.2) should be rewritten as 


UG Ug a a aa) i (= Ura); (7.3) 


and (7.3) now defines an approximation U7; to uj’, within the bounds of the dis- 
cretization and other errors. Note further that the representation of the numerical 
solution within the FDM framework is pointwise, i.e. we know pointwise values of 
approximations of the exact solution. 
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We conclude with a philosophical remark: 





In FDMs one discretizes the operators occuring in the differential equation. The 


discrete solution U” := (ur, li jgE€Z;ne N) consists of point values at the 


grid points of G” x G; j. 





In the FEM one always starts from a weak (variational) formulation like (2.2). The 
space Q := R? x {t € R | t > 0} is divided by means of a grid into certain elements. 
These elements may be chosen as the space-time bricks building the finite difference 
grid G” x G; j but can be far more general. Due to their simplicity and their ability 
to resolve even complex geometries, triangulations are often used where the elements 
are triangles or tetrahedra. Applying the notion of weak solutions and choosing a 
space W of test functions leads to 


V@ewWw: [ (uae + urð + udy®} dtdxdy=0. 


Replacing the test function space by a finite dimensional subspace W” leads to 


vo" ew: | {wae + ud," + ud,0") dtdedy =0. 
Q 


This is a set of so many equations as the dimension of W” := span{®,,...,®,}. A 
popular choice is the space of polynomials linear in t and x, y vanishing outside each 
element X C Q. 


Choosing the approximation 


U(a,y,t =5 ua, (x,y,t 
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for u then leads to a linear system for the unknown coefficients U; since 


Wg Shoe. yn: pO dt dx = 


[{Su0, a0; + 208.0, Dj DDT ®; 7 dt dx dy 
Q 


= i=1 i=1 


n 


= SoU; [08 dt dx dy + SOU; [vide dt dx dy 
Q E Q 


i=1 


+% U; [eae dtxdy | = 
i=1 Q 
=. J 08; + 260; + 0,05) at de dy U0. 
i=l Q 


Now, given U; on the old time level, the U; on the new time level follow from solving 
this linear system. 


FEMs of the type just described are usually called space-time FEMS. It is often conve- 
nient to discretize in space only with finite elements and tackle the time derivative in 
a finite difference manner. In order ro remove the coupling of every element with any 
other element one can use discontinuous Galerkin methods, where discontinuities 
are allowed at element faces. These types of methods are very close to finite volume 
methods, but not identical. Also popular are streamline diffusion FEMs which belong 
to the class of Petrov-Galerkin schemes. Petrov-Galerkin methods are characterized 
by choosing a test space different from the Ansatz space in which the discrete solution 
is sought. Although some mathematical analysis is available for streamline diffusion 
methods the numerical results obtained so far can not cope with results obtained by 
FDMs. 


Again we end up with a philosophical remark: 





In FEMs one discretizes the solution of the differential equation. The discrete 
solution is made up of polynomials defined piecewise on the elements. 





The FVMs as third class of methods share all good properties of the FDM as well as 
of FEMs. As with FEMs arbitrary cells in space can be used for the discretization, 
while all of the FD technology is applicable at the same time. Here, one starts with a 
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type of weak solution which stems directly from gas dynamics, where the conservative 
equations are formulated as integral balances on so-called control volumes. For our 
generic model this weak form results from integrating 


iu + ðzu + Oyu = jut+(V,u)=0, u:=(u,u)’, 


over a bounded control volume o C RÊ with boundary 0c where the unit outer nor- 
mal vector n is defined almost everywhere. From Gauss’ divergence theorem it then 
follows 


d 
£ | udedy+ fan) ds =0. 


Oo 


Now the cell average of u on ø is defined as 
= 1 
T(t) = Jo] feen t) dx dy ’ 
Oo 


where |a| denotes the area of ø, so that the weak formulation of gas dynamics enrolls 
itself as being an evolution equation for cell averages 


d 1 

LT, = —-— pe ds. 7.4 

ri ia] fu n) ds (7.4) 
Oo 


FVMs are derived from (7.4) by approximating the line integrals by a suitable Gauss 
formula. The product (u, n) is replaced by a so-called numerical flux function which 
is essentially a finite difference formula. 


Again we conclude with a philosophical remark: 


FVMs are discrete evolution equations for cell averages. There is freedom in 


representing the numerical solution. 





Some people like to view FVMs as being part of the Petrov-Galerkin family of FEMs 
with piecewise constant trial functions. The central idea in FVMs is to adapt FDMs 
to unstructured grids like triangulations and to use the computational advantages 
from FDMs and FEMs. 


There are even more discretization methods available than the three classes described 
here, but they may often be subsummed as particular instances of FDMs, FEMs and 
FVMs. A popular method in the computation of incompressible viscous flows, and 
now becoming more and more attractive for compressible fluid flow, is the class of 
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spectral discretizations. Here the Ansatz space consists of globally defined trigono- 
metric or polynomial functions on the whole domain (or, as in the spectral element 
method, on macro-elements). Choosing a FEM-like Ansatz leads to spectral Galer- 
kin methods. If the differential equation has to be satisfied pointwise at so-called 
collocation points one talks about spectral collocation or pseudospectral methods. 
These methods are very close to FDMs with large stencils. 


7.2 The Finite Difference Calculus 


Let us now study more intensively the structure of finite difference methods. 


Even one century before the invention of calculus through Leibniz and Newton, the 
foundations were laid for a discrete calculus employing differences of functions. 
Napier’s !!° and Briggs’ !!! invention of logarithms revolutionized naval navigation 
and necessitated interpolation formulae for computing intermediate values in tabu- 
lated data. It was the genius Thomas Harriot !!? who founded the field of finite dif- 
ference calculus. This field flourished in the days of Leibniz and Euler and came 
to enormous influence in the area of the numerical treatment of ordinary and partial 
differential equations. 


To give the reader again an idea of the spirit of the finite difference calculus, we 
consider the forward difference operator A defined by 


Au(2) := u(x + Ax) — u(x) , 


where Ax denotes an increment different from zero. It is often convenient to introduce 
the shift operator S through 


Su(x) := u(x + Az) 


so that A = S — I, where J is the identity operator. Assuming a smooth function 
x — u(x), Taylor series expansion results in 


u(x + Az) -53 — Axr” D”u(z), 


d d d 
where D” is defined as — o — o - - - o — and D? = I. 
dx dz dx 


V— times 
10John Napier (1559 - 1617); Merchiston 


"Henry Briggs (1556 - 1630); London, Oxford 
112Thomas Harriot (1560 - 1621); Oxford 
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Utilizing the shift operator, this is nothing but the representation 


Sulz) = (È Zaen) u(x), 
v=0 


where now, formally, the term in parathensis is exactly e4”, i.e. we have found the 
representation 

S = e^®P (7.5) 
and hence 

A= eS?) _ J. (7.6) 


Introducing the central difference operator ô by means of 
u(x) := u(x + Ax/2) — u(x — Az/2) , 


a further Taylor series 








1 1 
u(a + Ax/2) = 5° (+A) 5, D’u(a) 
v=0 ` 
reveals 
= 1 Ax Skal 

sul) = 2 T (SP) we), 
1.€. 

6 = 2sinh (=>) : (7.7) 


Inverting this equation gives a relation between the differential operator D with the 
central difference operator 6, namely 


2 ô 1 1? 1?3? 
D = sinh! = — | 8- — 8? + —— 8 —... 7. 
Ar 2 Az (s z3 t z5 ) 78) 


where 6” := 6060---06,. Equation (7.8) is an exact difference formula for the 


n—times 
operator D. 
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More equations like (7.8) and (7.6) can be constructed. 


We are mainly interested in the applicability of formulae like (7.8). To get an idea, we 
consider the simple linear transport equation pu + zu = 0 which we rewrite as 


Ou = L(D)u := —O,u, 


where the differential operator D is now partial derivation with respect to x. From 
(7.5), we deduce 


u(a,t-+ At) = (ae) ulz, t) 
and, since 0, = L( D), we conclude 


api (eA?) a (7.9) 


We now use (7.8) to eliminate D in terms of 6 in (7.9). This gives 


upt = (Comes eae) uP. (7.10) 


Note that (7.10) is not an approximate difference equation but an exact one. Retaining 
only the linear terms in the exponential function gives 


Uptt = (I + AtL(2/Azsinh ~1(5/2))) UP. 


Retaining also only linear terms in the expansion (7.8) and noting that L(D) = -D = 
—0,, leads finally to 


ô 

n+1 n 
U!? — IT — At— U 
i ( A ) i 


Expanding the central difference operator gives a finite difference method 


At 
Un 1 Ur n Ur 
i a i Ar ( i+1/2 ae) (7.11) 


for the transport equation. 


Hence, the finite difference calculus may be directly used to construct finite difference 
methods for partial differential equations. More constructions of this type can be 
found in the literature. 
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Now let the above difference scheme act on the true solution u. Employing again 
Taylor series expansions 


upt! = u(x,t) + Atõu(zx, t) + O(At?) 
n Az 1 Ax? 

Uni = ulz, t) + -z rulg, t) + aor Orla, t) + O(Az?) 
n Ax 1 Ax? 

Ue yj = ulz, t) - -z zula, t) + Hoy zule, t) + O(Az?) 





and inserting them in (7.11) gives 
2 At 3 
ut AtOju + O(At*) =u— a (Azd,u + O(Az*)) 


which, after rearranging, leads to 


ðu + Opu = O(At, Ax?) . (7.12) 


Thus, our simple difference method (7.11) is of order one in time and of order two 
in space. This notion of order is usually called truncation error since it results from 
truncating Taylor series. 


Note that we did not say a word about the convergence of (7.11). The property that 
(7.12) formally tends to ĝu + zu = 0 for At, Ax — 0 is called consistency. 


If a difference operator is applied to a sequence (U? | i € Z) then it acts on the se- 
quence index, i.e. A;U} := AU; = Uj", — UP. 


For later reference we introduce the antidifference operator A~! as follows. If 
Au; = gi holds, then we define 


u; = Ag +c, 


where c denotes an arbitrary constant. The constant appears since after application of 
A we get back Au; = g;. We first prove 


n-1 
X Aui = tn — ting » (7.13) 


i=Nno 


n—-1 
A (£ u) = Un. (7.14) 


i=n0 
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Equation (7.13) simply follows from 


n-1 
oy Au; = (tino+1 —Ung) +(Unp+2—Uno +1) +- . +(Un—1—-Un—2) +(Un—Un_1), 


i=no 
while (7.14) is nothing but 


n-1 n-1 n-1 n-1 n-1 
a(Sou)= E u-u Som Ym. 


i=no i+1l=no 1=NoO i=no—1 i=no 


Using formula (7.14), we can immediately deduce 
n-1 
An = So ute, (7.15) 
1=N0 
where no is an arbitrary start index. 


We now give a discrete analogue of the product rule, namely 


A(u;vi) = Su;Av; + viAu; . (7.16) 


The proof consists in evaluating the right and left hand side to see that they coincide. 


Writing the product rule in the form Au; = A(ujvi) — ui+1 Av; and applying the 
operator A~! results in A~!(vjAu;) = A7~'A(ujv;) — A7!(uj41Av;). Following 
formula (7.15), we finally get the formula for summation by parts, which is a discrete 
analogue to the formula of integration by parts, namely 


n-1 n-1 
X uAus = unon — So up Aute. (7.17) 
1=N0 i=No 


7.3 The CFL Condition 


In 1928 Courant!!? , Friedrichs''* and Hans Lewy published a milestone paper (cf. 
Section 6.1) which was forgotten in war times and became famous when it was redis- 
covered during the Manhattan project. 


'3Richard Courant (1888-1972); Göttingen, New York 
4K urt Otto Friedrichs (1901-1982); Brunswick, New York 
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As already emphasized in the last chapter, Courant’s interests at that time mainly con- 
cerned the proof of existence of solutions to partial differential equations. In that field 
he was one of the founders of the idea to employ finite difference schemes as approx- 
imations to partial differential equations. If he could prove that the difference scheme 
is consistent and convergent, and furthermore that the limit solution (as At and Ax 
tend to zero) satisfies the original differential equation, then he would have proven ex- 
istence of a solution of the partial differential equation. Although Courant, Friedrichs 
and Lewy succeeded with this scheme for the (elliptic) Laplace equation, they encoun- 
tered serious problems with the (hyperbolic) wave equation. Finally they succeeded 
in formulating the celebrated CFL condition which is necessary for convergence. 


In order to demonstrate the basic behaviour in the context under consideration, let us 
consider the simple forward expicit difference scheme for the slightly more general 
transport equation ĝu + aðru = 0, a < 0, namely 


At At At 
Use = U; = on ( i41 — Uj’) = (1 + ox) U; = an lin . (7.18) 





(7.18) is an explicit three-point scheme for the explicit computation of UIR, Viewed 
in the (x, t)-plane we can identify the difference stencil of (7.18) as shown in Fig. 
30. 


t 4 






































Figure 30: Difference stencil and characteristics 


The stencil gives rise to the definition of numerical characteristics which are straight 
lines enveloping the stencil. The numerical characteristics cut an interval off the axis 
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t = 0 which is called numerical domain of dependence. Changing the initial datum 
in this interval changes the numerical solution at point P while perturbations outside 
this interval would not affect the numerical solution at P. However, there is also the 
characteristic I of the differential equation we would like to approximate. If I lies 
outside the numerical characteristics as shown in Fig. 30, then changing the initial 
datum between IĮ and the numerical characteristic would change the true solution at 
P but not the numerical solution. Thus, convergence of the difference scheme is not 
possible. Keeping IĮ inside the numerical characteristics ensures that all perturbations 
of the initial datum which affect the true solution at P will always affect the numerical 
solution at P, too. 


The characteristic I has slope 


dt 1 
dra 
and the right numerical characteristic has slope —At/Az. Thus, necessary for con- 
vergence is At/Az < 1/|a| which is already the Courant-Friedrichs-Lewy (CFL) 


condition 
Az 


la| ` 


At < 


(7.19) 


In case of a > 0, and replacing the forward difference equation by the backward 
one, (7.19) again ensures the suitable behaviour of the domains of dependence. This 
dependence of the type of the difference equation on the sign of a, i.e. on the direction 
of the movement of the wave, leads to so-called upwind schemes. 


In the case of a scalar conservation law pu + ôx f(u) = 0, the same arguments hold 
for the quasilinear form ju + f’(u)0,u = 0. Hence, instead of (7.19), explicit finite 
difference schemes with three-point stencil have to satisfy the CFL condition 

At 


Ay mex |F (UL) <1 (7.20) 


in each time step n. 


7.4 Lax-Richtmyer Theory 


We now formalize the treatment of finite difference schemes supplementing chapter 6 
in order to present at least the historical central idea of a convergence theory for this 
particular type of schemes. 
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Although we have to confine ourselves with linear problems, the Lax-Richtmyer the- 
ory allows to gain deep insight into the interplay between the notions of consistency, 
stability and convergence. Originally, it was formulated for autonomous, linear initial 
value problems only but was later generalized by several authors so that also non- 
autonomous and nonlinear problems could be included!!». 


We consider an abstract evolution equation in a Banach space B, where a linear 
operator A : D4 — B acts upon a subspace D4 of B. The evolution equation we 
want to discretize is then given as 


Lut) = Aul), t€ [0,7], T> (7.21) 


u(0) = uo, ug € Dg. 


A one-parameter family u(t) with t € [0,7], u(t) € D4, is called a genuine solution 
of (7.21), if 














im etA- ut _ Au(t)]] =0. (1.22) 
At—0 At B 
O<t<T 


Let D be the set of all initial data up € B so that a unique genuine solution exists 
with u(0) = up and such that convergence in (7.22) is uniform in t!!6. Then there is 
a one-parameter family of linear operators Eo(t) : D — B defined by 


u(t) = Eo(t)uo 


which are called evolution operators or solution operators. We always assume that 
our problem is properly posed, which in particular implies that Fo(t) is bounded 
uniformly in t, i.e. there is a constant K > 0 sucht that ||Eo(t)||p 8 < K for 
0 < t < T. Then the Hahn-Banach extension theorem implies the existence of a 
generalized solution operator E(t) : B — B with the same bound K as F(t). The 
generalized solution operator has the semigroup property 


Vs,t>0: E(s+t) = E(s) E(t) 


and the closure of the operator A which we also denote by A is the infinitesimal 
generator of this semigroup. One can further show that E(t) commutes with A, so 





"Sef, R. Ansorge: Survey of equivalence theorems in the theory of difference approximations for 


partial initial value problems. In: Topics in Numerical Analysis (J.J. Miller, ed.), London-New York-San 
Francisco: Academic Press 1977 
116We shall see in a minute that uniform convergence is automatic. 
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that 
E(t) Aug = AE(t)uo. 


Then (7.21) can be written as 














E(At)-1 
lim po (FRO a) -0, 
At—0 At B 
0<t<T 


which, since || E(t)||B < K for0 < t < T, implies the uniform convergence of (7.22) 
with respect to t. 


In order to discuss fairly general finite difference schemes, let us introduce linear 
difference operators 


Most = 40/1 (At, Azi,..., Axa) 
and the class of difference schemes 
MUH = MoU”. 


We assume the existence of AS and that ATAO is a bounded linear operator defined 
in all of B. From the CFL condition we know that At and the Ax; can not go to 
zero independently of each other as far as explicit methods are concerned. Hence we 
assume relations Ax; = y;(At),i = 1,...,d, and introduce the discrete evolution 
operator 


C(At) = i (At, yi(At), ae) ya(At)) ` Ao (At, yi(At), ees) palAt)) 
so that the class of difference schemes under consideration is 


Unt! = C(At)U” . (7.23) 


The family of operators C(At) is said to define a consistent approximation for the 
initial value problem if for every u(t) in a class of genuine solutions with initial data 
u(0) dense in B the consistency condition 


(= 4) ue 


holds. Since u(t) denotes a genuine solution we can use the fact that due to (7.22) 

Au(t) is close to (u(t + At) — u(t))/At to give the alternative definition 

u(t + At) — C(At)u(t) 
At 


lim 


si (1.24) 
At—00<t<T 














B 


lim 
At—00<t<T 














=0. (7.25) 
B 
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The quantity 


u(t + At) — C(At)u(t) 


At = E(At, ug) 


B 














is called truncation error, and if this error is of order p, i.e. 
E(At, uo) = O(AP?) , 

the method is called a method of order p. Consistency therefore means 
E(At, ug) = o(1) 

which is particularly fulfilled for p > 0. 


Note that this abstract notion is in complete agreement with our earlier calculations in 
(7.12). 


Consistency guarantees that in the limit At — 0 the difference scheme formally re- 
produces the differential equation. The notion of convergence is concerned with the 
limit behaviour of the discrete solutions. Consistency does not necessarily already 
imply convergence! 


Due to (7.23), iterated application of C(At) gives 
U” = C(At)U" t} = C(At)oC(AthU™ 2 = --- = C(At)"U® = C(At)"uo , 


and this discrete function is expected to be a good approximation to u(nAt) = 
E(nAt)uo. The family of operators C(At) provides a convergent approximation 
to the initial value problem if, for fixed t € [0, T], 


lim ||C(Ajt)"4uo — E(t)uollg = 0 (7.26) 
Joo 


holds for each sequence (A;t) j=1,2,... tending to zero where n; is chosen to satisfy 
limjoo NjAjt =t. 


In reaching the discrete solution at time t = nAt, the operator C(At) has to act n 
times on the initial datum. Thinking of small errors in uo or errors due to rounding 
during the computation, we would not be satisfied with a difference scheme which 
would amplify these small perturbations until they spoil the solution. The central idea 
of stability is therefore to put bounds on the possible amplification of initial data. The 
operator C (At) is said to define a stable approximation if for some t* > 0 the infinite 
set of operators 


COCA, De aie?" O<SnAt<T, 
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are uniformly bounded, i.e. if there exists a constant « with 


|C(At)"Ip< k, VYneEN,VAt 
with 0< At< t andwith 0 <nAt<T. 


With the notions of consistency, stability and convergence we can now formulate the 
central Equivalence Theorem of the Lax-Richtmyer theory: 


Theorem 7.1 : A finite difference scheme approximating a properly posed initial value 
problem is convergent if and only if it is consistent and stable. 


PROOF (of the direction: consistency A stability > convergence!""): 
U(0) = uo implies 
U"—ult) = $ C(A) { C(At) u (tav) — u (tn—v41) } + (tn) u(t) 
v=1 
(ty = wAt). 


Thus, from the stability property, we obtain 
n 
|" -ule < Do CAD] - ICAL) un») = u (tav) 
v=1 


+ |lu(tn) — ule 


IA 


R Ss" ||C(At) Eo (tr—v) uo — Eo (tr—v+1) voll g 


v=1 
+ |lu(tn) = ule 
leading to 


JU" = utl < KnAt-é(At,uo) + |lulin) -= ulle - 





'7From an engineering point of view, this is the more important direction. The complete proof can 
be found in the original paper of Lax and Richtmyer cited in footnote 92. If consistency and stability 
are defined properly, theorems of the Lax-Richtmyer type can be proven even in the case of nonlinear 
difference schemes for nonlinear equations, as already mentioned in footnote 115. 


7.4 Lax-Richtmyer Theory 167 


Taking n At < T into account, we find 


JU" -ule < KT e(At,uo) + |lu(n At) —u)|lp - 


Because of n At — t , because of the continuity of u with respect to t, and because 
of the consistency condition (7.25), the truth of the assertion is obvious. 














Particularly forn At = t, we get an estimate for the global discretization error 
|U” — u(t)|| g namely 


JU" -ull <  Kte(At,uo) 


Here, round off errors are not included! 

Example: 

Let us again treat the transport equation ru + a0,u = 0 in the case of a < 0. 

If then (7.18) will be expressed in terms of the Lax-Richtmyer theory, and if a relation 
between the step sizes At and Ax will be described, we find 

At 
Ag ` 


At 
CAD < |i + aS] + la 


AT 
The CFL condition (7.19) then ensures 





IAD < (14 &) +E = 1, 


hence 


IC™ Allo < CAD S 1. 
Thus, the CFL condition guarantees the numerical stability of the scheme in the sense 
of the Lax-Richtmyer theory. 


And also the other direction is true: Stability in the sense of Lax-Richtmyer requires 
the CFL condition to be fulfilled, and this holds independently of the norm. 


Of course, the upwind scheme for a > 0 can be treated in the same way. 


The upwinding idea can also be looked upon from the point of view of information 
theory: Looking into the direction of the traveling information it seems to be unrea- 
sonable to use for the computation of a future value of the solution at a given space 
position such values of the present time level whose positions lie downstream and 
were not yet influenced by the information. 
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7.5 The von Neumann Stability Criterion 


Consistency of a finite difference scheme is easy to verify by means of Taylor expan- 
sions. In contrast, verifying stability using the pure definition may be much harder. 
In the case of linear initial value problems or problems with periodic boundary data 
there is a simple mechanism due to von Neumann!!8. 


We follow start with a complex Fourier transform of a discrete function represented 
as a sequence U” = (UP) jez and defined at grid points only. The discrete Fourier 
transform is defined as 


U Ey Ure, 0<€ <2r. 
jEZ 


The importance of the discrete Fourier transform in analyzing stability of finite differ- 
ence schemes can be seen if it is applied to the shift operator S defined already in sec- 
tion 7.2. For our purposes, we refine the definition by calling S}U” := S(Ax)U” := 
(U®.,)jez the forward shift operator. Then it follows from the definition of the 


j 
discrete Fourier transform 


TH) = Yup = Supe io 


jez jEzZ 
= P Upee = dY upei 
jez jEZ 
= FTE). 


Thus, applying the shift operator to a grid function is, in Fourier space, multiplication 
with e'S. The function ef$ is called the symbol of the shift operator 9. 


Completely analogously one gets for the backward shift operator S_U" 
S(—Azx)U” the symbol e~*, ice. 


mm 


S_UR(E) = e80 (e). 


Since every linear finite difference scheme is made up of combinations of forward 
and backward shift operators the way is open to look at discrete Fourier transforms of 
difference schemes. 





118John von Neumann (1903 - 1957) Berlin, Hamburg, Göttingen, Princeton 
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Example: 


Consider again the difference scheme 


U; ie U; = A i+1 7 Uj.) 


as an approximation to 0;u+a0,u = 0. It is an easy exercise to see that the scheme is 
consistent of order 2 in space and order I in time. Written in terms of shift operators 
and sequences we get 


aS (1 egress s.)) UF. 
x 
Now, taking Fourier transform on both sides, it follows that 


(E) = (1 = Aa (cit - E PE). 





If the symbol O(At, €) := 1 — Ata ( elf — e-i€) turns out to be larger than one in 
2AT 


modulus, this means that any small perturbation in U” will be amplified in the step 
computing U"*". In our case, since e+§ = cos € + isin €, we get 











At... 
O(At, €) =1— Are! sin€ . 


Taking the modulus, this amounts to 
ae 
|O(At, €)|? =1+ (Za) sin? é > 1 
Ax 


regardless how small the time step is chosen. The symbol is usually called the am- 
plification factor of the difference scheme. The present scheme is unconditionally 
unstable, since the modulus of the amplification factor is larger than one except for 
E € {0,7, 27}. Assuming a > 0, the upwind difference scheme 


At 
Cos ae (U; —UP1) 
results in the sequence formulation 
At 
unt} = | I — —a(I- S) ) U” 


which gives the discrete Fourier transform 


— 


Fg) = (1 - Fatt - 9) FE). 
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Let us define K := Sta and compute the amplification factor 
Q(At,é) = 1— K(1 =e) < |1 — K| + K le 4 
~w 
=1 
CFL condition 
1-K+Kk=1. 


Hence, our upwind scheme is stable under the usual CFL condition. 


We remark in passing that the same analysis holds in multiple space dimensions. Here, 
the discrete Fourier transform is defined as 


on =) => Ure ikTE= 
kezi 


where d is the space dimension and E = (&,... , Ea)T. 


In the case of systems of equations with n components in d space dimensions, ampli- 
fication is given by 


U+ (8) = @(At, 2) (2) 
where now © is an n x n matrix, the amplification matrix. For reasons of generality 
we stay with the case of systems in the following. 
Recursion reveals the relation 

U"(S) = 0” (At, 2)U(E) 


so that stability has to do with the properties of the iterated amplification matrix. In 
general, we can state the necessary and sufficient stability criterion that for some 
positive At* the matrices (or, in the case of scalar equations, the amplification factor) 


should be uniformly bounded for 0 < At < At*,0 < nAt < T andall E € Ze, 
If pọ(At, Z) denotes the spectral radius of © (i.e. the maximum of the moduli of the 
eigenvalues), then it follows from linear algebra that 


po (At, 3) < ]O"(At, &)|| < ]O(At,=) ||". 


A necessary condition for stability is therefore the existence of a constant C so that 


po(At,B) <C 
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for 0 < At < At*,0 < nAt < T and all & € Z°. Without loss of generality we 
assume C > 1 and therefore 


po(At,B) < Cr, O<n 


IA 


Thus, in particular, we get 
po(At,B) < CF. 


Now the expression C'^*/T is bounded for 0 < At < At* by a linear function 1 + 
const At. 


Thus, a necessary condition for stability is the celebrated von Neumann condition 
|Ai] < 1+ O(At) 


for0 < At < At*, al E € Z andi = 1,...,n, where A; denotes the i-th eigenvalue 
of O(At, B). 


It should be noted that there are several other stability conditions like the famous 
Kreiss matrix theorem or the much less known criterion of Buchanan! !’. 


7.6 The Modified Equation 


Finally, a tool should briefly be explained which can often help to analyse finite dif- 
ference procedures. 


To every finite difference approximation of order O(At?, Ax“) of a given differential 
equation there is another differential equation to which the difference scheme is a 
better approximation. 


In order to demonstrate this remark, let us start again from the simple upwind scheme 


for the transport equation 
Out+ad,u=0, a>O0. (7.28) 


"For a thorough treatment see: Richtmyer, R.D. and K.W. Morton, Difference methods for initial- 
value problems, 2nd ed., New York, London, Sydney: Interscience Publishers 1967 
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Inserting the Taylor series expansions of a sufficiently smooth function u, i.e. 








At? At? 
utl = u+ Atout => ou + bu + O(At*) 
Az? Az? 

ue, = u-—Axcdzut = opu — = ozu +O(AzÎ), 

into the upwind scheme results in 
At At? A Ax? 
Ou + > ou + Oe + O(At?) = —ad,u+ a= Ozu — a—— du 
+O(Az?) , 
or, after rearrangement of terms, in 
At A At? Ag? 
ðu +aðsu = -50u + D - =—dfu — = ~— u 


(7.29) 
+0(At?, Az?) . 


We do now want to replace the time derivatives on the right hand side by space deriva- 
tives. In order to do so, we are not allowed to use the differential equation (7.28), 
since our u does not necessarily satisfy this equation but equation (7.29)! Therefore, 
we have to accomplish our task using (7.29) alone. We take the time derivative of 
(7.29), i.e. 


aAx 
2 





At A 
Put aru = > Ou + 0,02u — su 


A 2 
2 = 3u + O(A, Az’) , 


and then the space derivative of (7.29) which we multiply by a, 





+ 


a Ar aAt? 


2 








At 
a0xu+ au = a= df On + 


aĉ Ar? 
6 


3u — oazu 





fu + O(At?, Az’) . 


Subtracting the last equation from the first one gives 





De DAO” Lag. Cis Gran GO ag 
Ou = a Ozu+ At ak 7 Ove + Ax gu 7 Ont 


+0(At?, Az?) . 
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Some further similar manipulations of (7.29) lead to 


Bu = -aðu + O(At, Az), 
0,02u = au + O(At, Az), 
0,02u = —ad?u + O(At, Az). 


Substituting these expressions into (7.29) yields the modified equation 








A At Az? (AP At 
Ou+adzu = a 1 — — u- $ 2 2—~ — 3— + 1 ou 
2 Ax 6 x x 


+O(Az?, Ax? At, ArAt?, At?) . 


What can be seen from this equation ? First of all we see that the truncation error of 
our difference scheme applied to the transport equation (7.28) is O(At, Ax), since 
the first term on the right hand side of the modified equation is 


sAadeu — 5 Ata 


which goes to zero linearly with time and space increments. Moreover, our upwind 
scheme is a good approximation to the parabolic differential equation 


ðu + aðru = c(a, At)02u 


with positive (due to the CFL condition) diffusion coefficient. Since such equations 
have smooth solutions, we expect our upwind scheme to be well-behaved. 


Many further properties of the difference scheme can be deduced from the modified 
equation, e.g. numerical stability. 


7.7 Difference Schemes in Conservation Form 


It seems to be obvious that numerical procedures can only be expected to lead to good 
approximations of the solution of the original problem if as many properties known in 
advance as possible are imitated by the scheme. 


One of these reproductions, already mentioned in chapter 6, is the use of so-called 
methods of conservation form in order to treat conservation law equations numeri- 
cally. 
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This means that the system (2.1) of conservation laws, i.e. 


AV + Orf(V) = 3V +I f(V) 3V =0 (7.30) 
is ean by an explicit (2k + 1)-step finite difference scheme 
1 
n+1 n Las n on = 
A V9) + (hy 93-4) =0 Ga 
with a consistent numerical flux 
fae = g (Vi ay Viako 2 8 Vo Vas Viak) (7.32) 


and not —e.g.— by a method like 


1 
a —(vitt_ve) 4 or) say (Vi - Vja) = 0. 


An example of a method of type (7.31) is the 3-point Lax-Friedrichs Scheme!”° 


n i 1 s 
vos 5 (V3 PVT) aed Vi) VO 4139) 
with 


A 1 Az A h 
Ij+44 = EDINI (Vi — V5 7) T3 5 {fl eee + £(VF)} $ 


If the numerical vectors V? in (7.33) are replaced by an exact smooth solution 
V (x£; , tn) of (7.30), the left side of (7.33) differs from zero, but its Taylor expansion 
then leads to a truncation error of order O(At) in the sense of the Lax-Richtmyer 
theory provided that A = at = const is prescribed. 


Hence, (7.33) is a method of conservation form of order 1. It can be shown that the 
method is monotone in the meaning of Section 6.4, and, as a matter of fact, all mono- 
tone methods are only of order 1 as can be established by means of the ‘Modified 
Equation’ of Section 7.6. 


But, of course, there are also methods of higher order. One of them is the 3-point 
Lax-Wendroff Scheme!”!. Its numerical flux reads in case of scalar 1-dimensional 
conservation laws of type (2.1) as 


4 
2 


1 1 
ag = 5 {POD +H) — 5 hg? a] (134) 


12 Proc. Nat. Acad. Sci. USA 68 (1971), 1686-1688 
121 Comm. Pure Appl. Math. 13, 217-237 (1960) 
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with 


we, = fad FO), _ ae 
ged — B.A 

IFZ vha = vp Az 

By Taylor expansion, it can be verified that this is a method of order 2. Therefore, 
the scheme seems to be a useful one, at least as long as the smoothness assumptions 
necessary for the Taylor expansion are guaranteed in a neighborhood of (x; , tn). 


But if the exact weak solution shows a shock at a certain position, the numerical 
solution begins to oscillate close to this position though the scheme respects the con- 
servation character of the original problem. 

But what is not respected by the method is the TVD-property of the solutions of the 
original problem. 


Numerical solutions produced by means of TVD-methods do not begin to oscillate 
because this would lead to an increase of the total variation. 


Several attempts were made in the literature to overcome the problem of the low order 
of monotone methods and the demand for keeping the advantages of higher order 
methods, e.g. the use of greater step sizes without loss of accuracy in areas where the 
exact solution is smooth. 


One of these ideas, e.g. in case of scalar problems, consists in the construction of 
so-called flux limiter methods: 


Let gy be the numerical flux of a higher order method, e.g. of the Lax-Wendroff 
scheme, and let gz be the numerical flux of a low order procedure, e.g. a TVD- 
scheme. 


If GH (UF 44-4 yo Vyk) is abbreviated by gy(v";7) and analogously gz(v”; 7), 
let us create a new method by 


(V5) = grlu” j) + plo” j) {gal 5) — rvs 5)} - 29) 


Herewith, y —called the limiter— has to be chosen in such a way that y ~ 1 in 
areas where the exact entropy solution is expected to be smooth, and y œ~ 0 close to 
shocks. In other words, p has to be a sensor for sudden strong increase or decrease 
of the unknown exact solution expected to occur at positions where the numerical 
solution increases or decreases strongly with respect to space, e.g. 


ve — yn 
y= (i=) (7.36) 
Vj41 Yj 
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where y is a sufficiently smooth and bounded function with y(1) = 1 and y(0) = 0. 
Of course, difficulties can still make trouble, e.g. if an extremum of the solution lies 
between x; and xj; +1. 


Also other ideas were created in order to get over the problem of non-physical oscilla- 
tions as far as higher order methods are under consideration. One of them in the case 
of higher dimensions was the introduction of so-called ENO schemes (essentially 
non-oscillatory schemes)!*”. This idea consists of a piecewise polynomial reconstruc- 
tion from the piecewise constant values of the approximate solution of the last time 
step where the stencil for each cell is adaptively chosen in such a way that intensive 
oscillations near shocks do not occur. 


7.8 The Finite Volume Method on Unstructured Grids 


Let us finally briefly put the idea of finite volume schemes on unstructured grids in- 
troduced in section 7.1 into concrete form. For convenience, we restrict ourselves to 
the scalar 2-dimensional conservation law (cf. (1.10)) 


iu + Or filu) + Oyfo(u) = 0 in R?xR*, 
(7.37) 
u(z,y;0) = uo(a,y) in R?. 


A so-called unstructured grid of Ê C RÊ? consists of a set of polygons o; (i € 
I c N), e.g. triangles, each of which has the same number m of edges and vertices, 
so that 


OL) ae 


iel 


and two polygons are assumed not to intersect or only to intersect along a common 
edge or at a common vertex. These polygons are often called cells. The angles at the 
vertices have to be bounded uniformly away from zero, and we assume h = O(At) 
where h represents a measure for the maximal diameter of the cells. 


Unstructured grids can more easily be adapted to complicated geometries than struc- 
tured grids, e.g. to the gas flow around flying objects, and refinements of the net in 
regions where the exact solution is expected to vary more than in neighbouring areas 


'2 cf. Harten, A.: Multi-dimensional ENO schemes for general geometries. ICASE Report 91-76, 
Langley Research Center: Hampton 1991 
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can easily be realized. The development of grid generation techniques, therefore, is 
a very important task in this context. 


Let I; C I be the index set consisting of the numbers of the particular cells being 
neighbours of g; with a joint edge. 


With u” (x,y) := u(x, y; n At) (n € No) in the case of a given time step At and with 
an exact solution u assumed for a short time to be smooth, equation (7.37) leads to 


aes y) = u” (x, y) 


AG + Or filu” (x, y)) + Oy fo(u (z, y)) = O(At) . (7.38) 


Taking the divergence theorem into account, integration of this equation over the cell 
g; with edges Oo; p (k = 1,2,---m) results in 


y s (fu), me) ds = O(At) 


1 n n 
xf {u Hx y)— u (x,y)} doi + 
Ti kel; 


(7.39) 


where f(u) = (filu), fo(u))” is assumed to be sufficiently smooth and where n't- 
denotes the outer normal unit vector on 0o;, jee ee 


If we replace the functions u” (x, y) , u"+! (a, y) along the cell o; by constant approx- 
imate values v? , v?™! to be calculated, if the term O(At) is neglected, if |o;| denotes 
the area of g; and s; y the length of the edge ôo; x, (7.39) leads to 





At 
1 ~ k 
eh a v = loil 3 (PER) , n” ) Sik (7.40) 
kel 
(n = 0,1,2,---) where the values ©, represent approximations for u along the 
edges ðc; k, respectively, and where v? is chosen to be the mean value of the initial 
function on øg; : 
0 1 À 
y= lal uo(z,y)do;, , Wiel. 
Gil 
Ti 


Obviously, from a physical point of view, the term (f(0",), n**¥) si describes ap- 
proximately the flux through the edge 0o;,x. 





121 
cf. (7.4) 
'22Here, we assumed that each of the cells has neighbours along each of its edges. If this is not the 
case, namely along the edges being parts of the boundary of Q, one has to respect suitable boundary 
conditions. 
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The question arises, how the values ù, should be chosen. 


Because the conservation principle should also be reflected by the approximate values, 
we aim for 


(Fn) > nF) sik = — (F (ki) n") spi (7.41) 


(Ski = Si,k). 


In order to create a scheme easy to handle, we are going to replace v;",, and vj, by 
a certain mean value of v? and vý together with a replacement of the physical flux 
by a numerical flux g; (V? , vj) that fulfills the relation (7.41) exactly, i.e. 


Gk (UE, Ue) = — GhilVee» Ui) - (7.42) 
Moreover, 
gi p(w, w) = (f(w), nF) six (7.43) 


should be fulfilled making the finite volume method now constructed, namely 





gers =v) — Ty Ji klU; , vg) (7.44) 
ail kel; 


consistent with (7.39). Finally, we do now forget the assumption on the smoothness 
of the exact solution because the integration process allows to deal also with weak 
solutions. 


Convergence proofs for FVMs in higher dimensions and on unstructured grids, par- 
ticularly mathematical proofs for convergence to the entropy solutions, are often not 
yet available. But what can be done is, to formulate the methods in such a way that 
they lead to well-known and well working methods if analogously formulated for 
one-dimensional problems. 


A method which could be shown to behave very well in the 1-D situation was the 
Engquist-Osher flux splitting scheme!? with a numerical flux described by formula 
(6.24). 


If the numerical flux g; x in (7.43) is chosen accordingly to (6.24) as 


gix(v,®) = [fin(o) + FO) (7.45) 





IB of. section 6.4 
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with 
Îi r(w) = (f(w), n™*) 5i,k 
and with 
v 
+(v) = f max { fl p(w), 0} dw + fix(0) 
0 
(7.46) 
A ô A 
fir) = f min { f p(w), 0} dw , 
0 
it corresponds to the Engquist-Osher scheme in 1-D, and it fulfills the relations (7.42) 
and (7.43). 


After these brief introductions, we stop here, but it should be mentioned that there 
are lots of schemes available for the numerical treatment of conservation laws like the 
Euler equations or for the solution of Navier-Stokes problems, scalar or systems, one- 
dimensional or multi-dimensional. Most of them show advantages compared with 
other schemes but often certain disadvantages as well. Many additional informations 
can be found in some of the extensive monographs listed behind this page. 
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